# LinStabPy > Sparse linear stability and resolvent analysis of compressible viscous flows. Python toolkit (CLI + library) that builds the linearized compressible Navier-Stokes operator on 2D structured base flows and solves eigenvalue and resolvent problems through PETSc/SLEPc with MUMPS sparse-direct factorization. LinStabPy answers two questions about a steady flow field you supply: stability (does a small disturbance grow on its own, at what frequency, with what shape?) and resolvent (at a given frequency, which harmonic forcing produces the largest response?). Base flows are loaded from files, never computed here. Cases are defined in YAML; results land in contract-stable `.npz` archives. How to use this file: find the task below, follow the link, and read that page. Links point at human-readable HTML pages; each page's Markdown source lives in the repo under `wiki/`. Critical facts agents get wrong: the `linstabpy` CLI takes NO case arguments -- config and output travel by `LINSTABPY_CASE_CONFIG` and `LINSTABPY_CASE_OUTPUT_NPZ` env vars, and unknown CLI bare words and `--flags` warn on stderr but are otherwise ignored (`parse_known_args`, exit code unchanged), so a guessed command still runs the wrong case. Install is the hard part (one MPI stack, complex-scalar PETSc/SLEPc + MUMPS via `install.sh --diagnose` -> `install.sh` -> `install.sh --check`); nothing that solves will work until `--check` passes. On macOS always use `.venv/bin/python`, never `uv run`. Fixed conventions repo-wide: temporal eigenvalues `omega = eig(i*L0)` with `exp(-i*omega*t)` time dependence, unstable exactly when `Im(omega) > 0`; SBP orders are 2, 4, 6, 8 only; default shifted solve is `KSPPREONLY` + `PCLU` through MUMPS. ## Start here - [Getting started](https://linstabpy.vercel.app/getting-started/): Install, first runs (example 0 needs no data), the example family, your first own case. Read this before running anything. - [Boundary interfaces](https://linstabpy.vercel.app/boundaries/): The machine surface: env-var CLI, YAML case schema, Python lazy imports, file inputs, NPZ outputs. Required reading for agents driving runs. - [Core concepts](https://linstabpy.vercel.app/core-concepts/): Exact definitions of base flow, growth rate, gain, Chu norm, wavemaker, sponge. Open beside any science page. - [FAQ](https://linstabpy.vercel.app/faq/): Which analysis do I need, shift-branch confusion, dissipation on/off, sawtooth wiggles, cache-valid surprise. ## Reference - [Overview](https://linstabpy.vercel.app/overview/): What LinStabPy computes, scope boundaries, system context, technical architecture. - [Architecture](https://linstabpy.vercel.app/architecture/): Packages, how they connect, main processes. Follows `docs/architecture.md` ownership map. - [Workflow](https://linstabpy.vercel.app/workflow/): Steps of a run in order: config to baseflow to mesh to L0 + BCs to factorize to eigen/resolvent to payload. - [Module guide](https://linstabpy.vercel.app/modules/): Per-package quick reference plus deep package notes (driver, operators, mesh, solvers, router, linalg). - [Solver glossary](https://linstabpy.vercel.app/solver-glossary/): Every KSP, PC, EPS, ST, and MatShell term defined from zero with official PETSc/SLEPc links. ## Optional - [References](https://linstabpy.vercel.app/references/): Papers and manuals behind the wiki, grouped by what they justify. - [Repository README](https://github.com/openfluids/linstabpy#readme): Install matrix, MPI rules, example runner, troubleshooting. Same content as Getting started, repo-side. --- # Full content The pages linked above, inlined in order. Generated by tools/generate_llms_full.py — do not edit. --- ## Getting started Source: https://linstabpy.vercel.app/getting-started/ # Getting started From a fresh machine to your first spectrum and your first [gain](../core-concepts.md): install the stack, run the examples in the right order, read the outputs, then start your own case. Assumes linear algebra, [MPI](https://grokipedia.com/page/Message_Passing_Interface), and a Unix shell; no prior experience with this codebase needed. Time budget: install takes 30 to 90 minutes ([PETSc](https://en.wikipedia.org/wiki/PETSc) and [SLEPc](https://en.wikipedia.org/wiki/SLEPc) compile from source). Example 0 runs in under a minute. The [Blasius](https://grokipedia.com/page/Blasius_boundary_layer) example takes a few minutes on a workstation. The order that works: :material-numeric-1-circle: Install with `install.sh` (or conda if you already live there). :material-numeric-2-circle: Run example 0. It needs no external data and proves the numerics import chain. :material-numeric-3-circle: Fetch the LinStab2D reference data, run example 1, read its plots. :material-numeric-4-circle: Run the rest of the family through the runner script. :material-numeric-5-circle: Write your own YAML case and run it through the driver or a module. Math conventions used on this page: temporal [eigenvalues](https://grokipedia.com/page/Eigenvalues_and_eigenvectors) are $\omega = \mathrm{eig}(i L_0)$ with $\exp(-i\omega t)$ time dependence, so $\mathrm{Im}(\omega) > 0$ means growth. The [resolvent](https://grokipedia.com/page/Resolvent_formalism) is $H = C F (L_0 - (\gamma + i\omega) I)^{-1} F B$. You need these two lines to read any spectrum or gain plot below. ## What you need before you start **Operating system.** Linux is the reference path. macOS works with one [BLAS](https://en.wikipedia.org/wiki/BLAS) caveat (see [macOS and clusters in brief](#macos-and-clusters-in-brief)). Windows is not supported. The installer expects a POSIX shell with `bash`. **Compilers and tools.** C, C++, and Fortran compilers, plus `make`, `cmake`, `curl`, and `pkg-config`. One MPI stack with development headers (`mpicc`, `mpicxx`, `mpifort`, `mpiexec`). Optional: HDF5 development files, only if you want the MATLAB v7.3 parity tooling. **The one rule that saves hours.** One MPI stack per install. PETSc, SLEPc, `mpi4py`, `petsc4py`, `slepc4py`, and the `mpiexec` you run with must all come from the same MPI family. Mixing Open MPI headers with an MPICH launcher is the most common broken install. `install.sh --diagnose` exists to catch exactly this before anything compiles. **Disk and time.** The build directory (`.venv`) holds a full PETSc and SLEPc. Expect several GB and a long compile on first install. Later installs reuse the pattern and go faster. ## Install path A: the reference script `install.sh` is the supported workstation build: default Python 3.13, PETSc 3.24 with [MUMPS](https://en.wikipedia.org/wiki/MUMPS_(software)), complex scalars, editable install of this repo into `.venv`. Run four commands in this order. ```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 ``` What each step proves: | Step | What it does | Success looks like | |---|---|---| | `--diagnose` | Read-only host preflight. Checks compilers, MPI launchers, headers, dual-MPI conflicts, and that `mpiexec -n 1 /bin/true` works. Installs nothing. | :white_check_mark: No errors. If it complains, fix the host first. | | `install.sh` | Creates `.venv`, builds or installs Python with PETSc and SLEPc, repairs `mpi4py` when the wheel ABI mismatches your MPI, installs this repo editable. | :white_check_mark: Ends without error. It never calls `sudo` unless you opt in (see knobs below). | | `--check` | Verifies `petsc4py` and `slepc4py` import, PETSc uses complex scalars, and a tiny SLEPc non-Hermitian Krylov-Schur solve works from Python. | :white_check_mark: All probes pass. A real-scalar PETSc or a failed NHEP probe means the env is unusable for this code. | | `benchmark.sh` | Runs a small benchmark (example 0 plus a small example 1) and writes artifacts under `benchmark_runs//`. | :white_check_mark: A fresh run directory with a summary CSV and logs. | Useful variants: ```bash bash install.sh --force # wipe .venv and rebuild everything bash install.sh --force --latest # same, but newest compatible releases instead of pins ``` !!! warning "Read the log, do not guess" A healthy install prints which MPI launcher configured PETSc and SLEPc, which launcher verified them, the resolved MPI library directory, and single-rank and multi-rank check results. If those paths point at different MPI families, wipe and rebuild from a clean environment. Do not paper over a mixed MPI install with `LD_LIBRARY_PATH`. Knobs you may need (all in the README with full detail): | Variable | When to set it | |---|---| | `LINSTABPY_AUTO_INSTALL_SYS_DEPS=1` | Let the script call your package manager for missing Linux system tools. | | `LINSTABPY_SKIP_SYS_DEPS=1` | Skip system dependency checks on provisioned cluster nodes. | | `LINSTABPY_MPI_STACK=intelmpi` | You already loaded Intel MPI and want it used. | | `LINSTABPY_PETSC_MPIEXEC=/path/to/mpiexec` | PETSc configure hangs or picks the wrong launcher. | | `LINSTABPY_VERIFY_MPIEXEC=/path/to/mpiexec` | Install-time checks need a different launcher than configure. | | `LINSTABPY_MPI4PY_NO_BINARY=1` | Force `mpi4py` to build from source against your active MPI. | | `LINSTABPY_WITH_H5PY=1` | Also install `h5py` for MATLAB v7.3 file tooling. | Three launcher knobs, three jobs: `LINSTABPY_PETSC_MPIEXEC` picks the launcher PETSc/SLEPc configure probes, `LINSTABPY_VERIFY_MPIEXEC` picks the launcher install-time checks run under, and `LINSTABPY_MPIEXEC` picks the launcher `run_all_examples.sh` uses. Each defaults to the auto-detected working launcher; set them only when the probe picks the wrong MPI family. `benchmark.sh` re-exports its resolved launcher as `LINSTABPY_MPIEXEC` (and its interpreter as `LINSTABPY_PYTHON_BIN`) for its child runs. Run everything through `.venv/bin/python` directly (scripts should do this so logs show the interpreter). For interactive shells, `source .venv/bin/activate` is convenient. ## Install path B: conda Pick conda if you already manage environments that way and want to skip the PETSc compile. conda-forge supplies Python with a complex-scalar PETSc, SLEPc, and `mpi4py` stack. pip then installs only this package, editable, with `--no-deps` so it never pulls a second NumPy or compiles PETSc. ```bash # Unset these first. A leftover PETSc prefix binds petsc4py to the wrong build. unset PETSC_DIR PETSC_ARCH SLEPC_DIR conda env create -f environment.yml conda activate linstabpy python -c "import numpy as np; from petsc4py import PETSc; assert np.issubdtype(PETSc.ScalarType, np.complexfloating)" ``` That assert is the whole gate. If it fails, the environment cannot run this code even though `import linstabpy` still works (plain import never touches PETSc by design). Constraints to know: * Python is `>=3.10,<3.14` here. Python 3.14 stays an `install.sh` freshness test; conda-forge `petsc4py` may not resolve it. * Never add `petsc4py`, `slepc4py`, or `mpi4py` to the pip list. Those sdists compile PETSc from source and fight the conda MPI stack. * `environment.yml` is not a second reference build. If conda-forge PETSc differs in minor version from the `install.sh` build, record both versions with your results. * There is no `environment-gpu.yml`, deliberately. GPU PETSc is not a supported conda path. See [What is still missing](#what-is-still-missing) for the current GPU status. ## The CLI as it really is This section corrects a mistake that lived in older docs. The `linstabpy` command takes **no case arguments**. It reads everything from environment variables and caches results by config hash. From `src/linstabpy/driver/cli.py`, the only flag is `--version`. ```bash linstabpy --version # light: prints version, never initializes PETSc or MPI linstabpy # full run; config and output come from env (below) python -m linstabpy.driver.cli # same entrypoint as a module ``` | Variable | Default | Meaning | |---|---|---| | `LINSTABPY_CASE_CONFIG` | `./configs/parity_config.yaml` | YAML case file. Legacy aliases `LINSTABPY_CASE2_CONFIG`, `LINSTABPY_CASE1_CONFIG` also work. | | `LINSTABPY_CASE_OUTPUT_NPZ` | `{prefix}_{config_md5}_{git_sha8}.npz` in cwd | Where the result archive lands. Per-run aliases exist too. Set a unique path per concurrent run. | | `LINSTABPY_CASE_LOG_PREFIX` | `linstabpy` | Prefix for log lines and the default output stem. | | `LINSTABPY_FORCE_RECOMPUTE` | unset | Set to `1`, `true`, `yes`, or `on` (case-insensitive) to skip the cache check and recompute. CLI gate only; library `driver.run` callers always compute. | | `LINSTABPY_STAGE_HEARTBEAT_SECONDS` | `30` | Progress heartbeat interval for long stages. | A typical MPI run: ```bash export LINSTABPY_CASE_CONFIG=cases/my_case.yaml export LINSTABPY_CASE_OUTPUT_NPZ=results/my_case.npz mpirun -n 4 linstabpy ``` !!! warning "Unknown arguments warn, then are ignored" The CLI uses `parse_known_args`, so `linstabpy run case.yaml` does not error. It prints a warning on stderr (rank 0 only) — for bare words and `--style` options alike — then runs with the env config and drops your words on the floor. If a run used the wrong case, this is the first thing to check. Single-dash PETSc options (plus one following value, e.g. `-ksp_rtol 1e-8`) still pass through to the solver stack. Caching behavior: on rank 0 the driver loads the existing output NPZ if present and compares its stored config hash and git SHA against the current ones. Match means it prints `Cache valid` and skips the compute. Mismatch or unreadable file means it recomputes. In parameter sweeps, either point each run at its own output path or export `LINSTABPY_FORCE_RECOMPUTE=1` deliberately. Most examples in `examples/` do **not** go through this CLI. They run as Python modules with their own YAML files (next sections). The driver CLI is for packed-field cases and campaign use. ## macOS and clusters in brief **macOS.** Same four commands as Linux. One extra step first when possible: ```bash brew install openblas scalapack ``` Reason: PETSc-downloaded ScaLAPACK can resolve [LAPACK](https://en.wikipedia.org/wiki/LAPACK) symbols to Apple Accelerate inside a Python process, which aborts Python `slepc4py` NHEP solves inside Accelerate `ZLADIV`. The installer prefers Homebrew OpenBLAS and ScaLAPACK at configure time and wraps `.venv/bin/python` so the right BLAS loads first. `install.sh --check` runs a tiny NHEP probe that catches this before examples run. On macOS, always use `.venv/bin/python` directly, never `uv run`, which can bypass the patched layout. **Clusters.** Load compiler and MPI modules first so the scripts see the same `mpicc`, `mpicxx`, `mpifort`, and `mpiexec` you will run with: ```bash module load bash install.sh --diagnose bash install.sh bash install.sh --check ``` For Intel MPI with oneAPI, source `setvars.sh`, export `LINSTABPY_MPI_STACK=intelmpi` plus the `I_MPI_CC/CXX/F90/FC` compiler variables, then run the same sequence. Full commands live in the README install section. **Thread policy for every MPI run.** One thread per rank unless you are benchmarking deliberately: ```bash export OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 export NUMEXPR_NUM_THREADS=1 VECLIB_MAXIMUM_THREADS=1 ``` Without this, each rank spawns a full BLAS thread pool and [sparse](https://grokipedia.com/page/Sparse_matrix) direct solves get slower, not faster. More ranks are not automatically faster for factorization either. Compare rank counts on your machine. ## First run: example 0, no data needed Example 0 validates curvilinear derivatives. It is pure NumPy, runs single-process, and needs no reference files. Run it first. If it fails, the problem is your Python environment, not MPI or data. ```bash .venv/bin/python -m examples.martini.example_0_curvilinear_validation ``` Outputs, next to the script: `example_0_curvilinear_validation.npz` (numbers), `.log` (run log), and PNGs including the curvilinear mesh figure and a field-derivative triptych. Open the triptych. It shows the test field and its discrete derivatives on the deformed mesh. This is also your first look at the [SBP](https://grokipedia.com/page/Summation_by_parts) operators in action: the derivatives stay clean at boundaries, which is the whole point of the validated closures. Re-plot any time without recomputing: ```bash .venv/bin/python -m examples.martini.example_0_plots ``` ## Second run: Blasius resolvent Example 1 is the real first science run: resolvent gains with optimal forcing and [response modes](../core-concepts.md) for a [compressible](https://grokipedia.com/page/Compressible_flow) Blasius [boundary layer](https://grokipedia.com/page/Boundary_layer) ($Re = 6.0\times 10^5$, $Ma = 0.3$, SBP order 4, grid $300 \times 100$). It needs one external file, the MATLAB base flow `BL_BF.mat`, which this repo does not bundle. Fetch the reference data once: ```bash git clone https://github.com/LinStab2D/LinStab2D.git ~/linstab2d-reference export LINSTABPY_BASEFLOW_DATA_DIR=$HOME/linstab2d-reference ``` Without that variable, the example raises `FileNotFoundError` telling you to set it. Put the export in your shell profile if you run these often. The Song and Savarino cases ship their own data or fetch scripts and never read this variable. Run the case, then the plots: ```bash .venv/bin/python -m examples.martini.example_1_blasius .venv/bin/python -m examples.martini.example_1_plots ``` What you get: * `example_1_blasius.npz`: gains, forcing and response modes, mesh, [Chu](../core-concepts.md) weights, BC provenance, solver metadata. Load with `numpy.load`. The key names form a compatibility surface documented in `docs/runtime_npz_contract.md`. * `example_1_blasius_baseflow.png`: the mean flow you linearized about. Always look at this before trusting modes. * `example_1_blasius_modes_forcing_response_{u,v}[_imag]_mosaic.png`: one mosaic per component. Each shows the optimal forcing (what to inject), the response (what grows), and their imaginary parts. How to read the mosaic: the forcing panel is the input shape $v$, the response panel is $u = Hv/\sigma$, and the gain $\sigma$ printed on the figure is their energy ratio in the Chu norm. A large gain with a clean, smooth mode pair is the expected picture. Noise or grid-scale wiggles point at resolution, filter, or boundary setup, not at the solver. ![Blasius base flow: streamwise velocity with the boundary layer developing downstream](figures/example_1_blasius_baseflow.png) *The mean flow you linearize about. Always look at this before trusting modes.* ![Blasius optimal forcing (left) and response (right) for the streamwise velocity component](figures/example_1_blasius_modes_forcing_response_u_mosaic.png) *What a healthy answer looks like: smooth, cleanly separated forcing and response structures. Compare your run against this before tuning anything.* The case definition lives in `examples/martini/example_1_blasius.yaml`: physics, grid, sponge, frequency, mode count, BC recipe, solver policy, each with a comment. YAML keeps the explanation next to the knob, which is why the repo uses it instead of long flag lists. [Your first own case](#your-first-own-case) walks a trimmed version. ## The example family map Six Martini examples form the core family. Run one, run all, then compare against MATLAB: ```bash bash run_all_examples.sh --list # what exists bash run_all_examples.sh --example 1 --skip-parity # one example, no MATLAB check bash run_all_examples.sh # everything plus parity and plots ``` | # | Module | What it teaches | Needs | |---|---|---|---| | 0 | `example_0_curvilinear_validation` | SBP derivatives on deformed meshes | Nothing. Pure NumPy. | | 1 | `example_1_blasius` | Resolvent gains, forcing and response | `BL_BF.mat` via data dir | | 2 | `example_2_jet` | Resolvent on a turbulent jet base flow | `M09jet.mat` via data dir | | 3 | `example_3_streaky_bl` | Temporal and spatial EVP plus Floquet options | MATLAB fixtures for plots | | 4 | `example_4_cavity` | Direct and adjoint eigenvalues, spectra | `Step_BF.mat` via data dir | | 5 | `example_5_paraboloid` | Resolvent on a 3D body-fitted mesh | `Paraboloid_BF.mat` via data dir | Each example has a matching `example_N_plots.py` that renders figures from the NPZ without recomputing. The runner calls compute then plots automatically. Direct module invocation writes the NPZ only. MPI form, with PETSc options passed after `--`: ```bash bash run_all_examples.sh --mpi -n 10 -m examples.martini.example_3_streaky_bl bash run_all_examples.sh --mpi -n 1 -m examples.martini.example_1_blasius -- -eps_monitor ``` Without `-n`, the runner defaults to half the detected physical cores. Override the launcher with `LINSTABPY_MPIEXEC` and the rank count with `LINSTABPY_MPI_RANKS` when the machine needs explicit paths. Interpreter: `LINSTABPY_PYTHON` wins over `LINSTABPY_PYTHON_BIN`, else `.venv/bin/python`. Beyond Martini, three case families cover harder ground. Leave them for week two: | Family | What it is | Data story | How to run | |---|---|---|---| | `examples/song/` | Song and Robinet SWBLI global stability, $\theta = 32°$ | Ships packed `song_32deg.npz` for checks; full pipeline needs Plot3D `grid.grd` + `mean_flow.dat` in `inputs/raw/` | `prepare_baseflow`, `plot_baseflow`, then `mpiexec -n 4 python -m examples.song.song --config ...` | | `examples/sawtooth/` | Lui 2022 suction-side SWBLI, stride-2 field | Ships packed `sawtooth.npz` ($300 \times 340$) | Baseflow figure, `mpiexec -n 4` EVP run, spectrum figure. See its README. | | `examples/savarino/` | 30.8° SWBLI resolvent gain maps and backend benchmarks | `fetch_zenodo.py` downloads the Tecplot base flow | Slurm campaign scripts. Validation and speed reproduction, not new physics. | | `examples/prescribed_forcing/` | Forcing demos | Local | Small demos, see directory. | Parity checks compare Python results against MATLAB reference exports: ```bash .venv/bin/python examples/martini/compare_martini_examples.py ``` Missing reference bundles skip cleanly here (silent in CI) and validate when present locally. The pytest parity marker behaves more strictly; see [Verify with tests and benchmarks](#verify-with-tests-and-benchmarks). ### What things cost No maintained per-example runtime table exists yet, but two measured anchors calibrate expectations. Both come from the Savarino case on one c9 node (dual Xeon Silver 4210R, 20 cores, MKL), August 2026, recorded in `examples/savarino/*.csv` with commit, ranks, and Slurm job per row. Single resolvent point, Mack point, 450,000 unknowns, 23.4M nonzeros, MUMPS (strong scaling, jobs run serially): | Ranks | Wall (s) | Peak RSS (GiB) | Gain | |---|---|---|---| | 4 | 190.5 | 18.6 | 179.1173862845 | | 8 | 154.6 | 25.4 | 179.1173862854 | | 16 | 148.2 | 40.6 | 179.1173862846 | Read it as: gains identical to 10 digits, ranks 4 to 16 buy 22 percent less time for twice the memory. Factorization dominates; more ranks are not automatically faster. Full 200-point gain map, same hardware, 5 jobs concurrent on 20 cores: 83.6 minutes as independent per-point jobs, 59.7 minutes through the line driver with anchor reuse. The router earns its keep on dense grids. Measure your own runs the same way: `benchmark.sh` writes wall times and logs under `benchmark_runs//`, and every result NPZ carries `cost_wall_*` keys plus provenance. When reporting a slow machine, send those plus the `--diagnose` output. ## The tutorial PDF `docs/tutorial/tutorial.pdf` walks the Martini examples figure by figure. Section 5 explains how to read the baseflow, spectrum, and mosaic figures. The PDF is not tracked in git, so build it: ```bash cd docs/tutorial && make tutorial.pdf # needs pdflatex ``` `make clean` removes build artifacts. If you lack a TeX distribution, the example PNGs plus [Second run: Blasius resolvent](#second-run-blasius-resolvent) cover the same ground without the prose. ## Your first own case You have two ways to run custom work. Example modules for MATLAB-parity style studies with generated meshes. Packed-field cases through the driver for your own base flows. This section sketches the packed-field route, which is the one new research takes. **A case is a YAML file plus baseflow arrays.** The YAML names physics, grid, sponge, boundary recipe, analysis targets, and solver policy. The baseflow is eight arrays with fixed names: `x`, `y`, `rho`, `u`, `v`, `w`, `T`, `sponge`. The loader reads either an `npz` archive or a directory of `npy` files. It never guesses filenames; the YAML names the file. Trimmed from the Blasius YAML, the shape of every case: ```yaml Re: 6.0e5 # Reynolds number, rho_ref U_ref L_ref / mu(T_0) Ma: 0.3 # Mach number, U_ref / sqrt(kappa R T_ref) Pr: 0.7 # Prandtl number kappa: 1.4 # heat capacity ratio (gamma in most textbooks) T_0: 273.15 # freestream static temperature [K] nx: 300 ny: 100 fd_order: 4 # SBP interior order: 2, 4, 6, or 8 sponge: amplitude: 5.0 xs_trans: 0.02 # ... ramp widths and offsets that keep outgoing waves from reflecting bc: penalty: 1.0e4 recipe: - {type: dirichlet, borders: lrbt, variables: uvwT} solver: # eigen targets, tolerances, backend selection ``` Three details that bite newcomers: 1. **Nondimensionalization is convective.** Density, velocity, and temperature equal 1 in the freestream. Pressure never appears explicitly; the nondimensional [ideal gas law](https://grokipedia.com/page/Ideal_gas_law) is $p^* = \rho^* T^* / (\kappa Ma^2)$. If your base flow comes from an acoustically scaled LES code, convert first. The README nondimensionalization section has the exact factors. 2. **Boundary recipes are data.** Each rule names a type (`dirichlet` or `neumann`), border letters (`l`, `r`, `t`, `b`, combined like `lt`), and variable blocks (`u`, `T`, `ruvw`, `ruvwT`). [Neumann](https://grokipedia.com/page/Neumann_boundary_condition) rules usually want `exclude_corners: true` and, on deformed meshes, `curvilinear_normals: true`. Assembly applies BC rows last and publishes the free degree sets downstream solvers must use. The `penalty` beside the recipe pins Dirichlet DOFs (diagonal value after row and column zeroing) and scales Neumann rows; keep the default unless leakage diagnostics complain. Full mechanics in [core concepts](../core-concepts.md). 3. **The dissipation epsilon has two halves.** Building the dissipation operator and setting its strength are separate flags, and setting only one is a silent no-op. The Workflow page documents the exact pair. On rectilinear meshes you usually want it off; on meshes curved in both directions you usually want it on. Run the case through the driver ([The CLI as it really is](#the-cli-as-it-really-is) for the env vars), then inspect the archive: ```python import numpy as np d = np.load("results/my_case.npz") print(d["omega"]) # temporal eigenvalues near your target print(d["gains"]) # resolvent gains, if a resolvent stage ran print(d["provenance_config_hash"], d["provenance_linstabpy_git_sha"]) ``` Keep the YAML next to the NPZ. That pair, plus the git SHA inside the archive, is what makes the run reproducible and is what the stability acceptance postchecks (`docs/runtime_npz_contract.md`) validate against. Before trusting any spectrum, check the basics the README lists in full: requested `nev` versus converged `nconv`, subspace size `ncv`, [EPS](../solver-glossary.md) tolerance and convergence test, postcheck residuals, boundary leakage, and target distances. Converged Ritz pairs with bad full-operator residuals happen on interior non-Hermitian problems, and the postcheck is what catches them. ## Verify with tests and benchmarks The test suite has three lanes: ```bash .venv/bin/python -m pytest --collect-only tests # what exists .venv/bin/python -m pytest tests/ # default lane LINSTABPY_MATLAB_PARITY=1 .venv/bin/python -m pytest -m matlab_parity # MATLAB comparisons ``` Deliberate semantics: :white_check_mark: parity tests are **deselected** by default, not skipped, so they never pad the skip count. Every skip in a default run names an absent optional backend ([STRUMPACK](https://portal.nersc.gov/project/sparse/strumpack/), [CUDA](https://en.wikipedia.org/wiki/CUDA), `superlu_dist`) on your machine. But once you request parity with the env var, a missing reference is a :x: **failure**, not a skip. A requested comparison that silently does nothing is worse than useless. Tracked files must stay under 60 MiB each: ```bash .venv/bin/python -m pytest tests/test_tracked_file_size.py -q ``` Large generated outputs, raw MATLAB exports, and machine-local benchmark products stay out of git. Regenerate or archive them externally. `benchmark.sh` artifacts live under `benchmark_runs//` with a summary CSV and logs. Send those plus full `--diagnose` and install output when reporting a broken machine. ## Troubleshooting quick map | Symptom | Likely cause | First check | |---|---|---| | MPI compile failure, `mpi.h` missing | Missing dev headers | Install your MPI stack dev package (for example `libopenmpi-dev`), or load the full cluster module. | | PETSc built against one MPI, launched with another | Dual-MPI environment | Pick one stack, put it first in `PATH`, wipe `.venv`, reinstall. Never mix. | | `mpi4py` import fails after switching MPI | Wheel ABI mismatch | Wipe `.venv`, then `LINSTABPY_MPI4PY_NO_BINARY=1 bash install.sh`. | | Python NHEP solve aborts on macOS, C SLEPc works | Accelerate LAPACK interposition | `brew install openblas scalapack`, reinstall, use `.venv/bin/python` only. | | PETSc configure hangs | Wrong launcher probed | `LINSTABPY_PETSC_MPIEXEC=/path/to/working/mpiexec bash install.sh`. | | Example raises `FileNotFoundError` for `.mat` | Data dir unset | Clone LinStab2D, export `LINSTABPY_BASEFLOW_DATA_DIR`. | | `linstabpy` ran the wrong case | Args ignored (stderr warning, exit 0) | The CLI takes no case args. Set `LINSTABPY_CASE_CONFIG` and rerun. | | Stale results reused | Cache hit | Check for the `Cache valid` line. Set unique output paths or `LINSTABPY_FORCE_RECOMPUTE=1`. | | Factor crash or OOM | Backend, ranks, or memory | Verify backend build, raise ranks, consider BLR options, check BC setup. | The README troubleshooting section carries the full commands for each row. ## What is still missing Honest gaps, so you can plan around them: * **No GPU install path.** `cudss` (through `nvmath-python`, wrapped around [FGMRES](https://grokipedia.com/page/Generalized_minimal_residual_method)) is the validated GPU route and matches MUMPS to 1e-12 on tested cases, but there is no one-command GPU environment. `install.sh --stack gpu-strumpack --dry-run` only prints a plan. CPU with MUMPS stays the default you should build on. * **No Windows support.** POSIX shell and Unix MPI are assumed throughout. * **Reference `.mat` files are external.** The four Martini base flows live in a LinStab2D clone you fetch yourself. Nothing in this repo vendors them, for provenance and size reasons. * **No per-example runtime table yet.** Wall times depend on grid, ranks, and backend. [The example family map](#the-example-family-map) has Savarino reference anchors (single point scaling plus a 200-point map) and the measurement protocol; Martini per-example timings are still unmeasured contributions welcome. * **Tutorial PDF needs a TeX toolchain.** `pdflatex` must be installed to build it. There is no hosted HTML version of the tutorial. * **Spatial and map entry points are library-first.** Temporal EVP and single-frequency resolvent have the smoothest CLI story. Spatial solves and $\beta$ and $\omega$ maps are best driven through module entry points and the `resolvent_map` API today. * **Error messages assume context.** Backend probing fails loudly (good), but some PETSc-level errors still surface raw. Paste them with the benchmark artifacts when asking for help. ## Where to go next * Run the family: [example family map](#the-example-family-map) commands, then the parity comparison. * Read the [Overview](../overview.md) for scope and vocabulary. * Read the [Workflow](../workflow.md) before writing your own YAML. It documents stage order, the dissipation pair, and acceptance gates. * Open one [module guide](../modules/index.md) package note when you need math or code detail. Start with LNSE Operator Assembly or Resolvent Analysis. * Keep the [Boundaries](../boundaries.md) page open while building cases. It lists every env var, flag, and output key. * Open the [solver glossary](../solver-glossary.md) beside any solver output. * Citing results: see [Cite linstabpy](../index.md#cite-linstabpy) and the [references](../references.md) page. The glossary defines `KSP`, `PC`, `EPS`, `ST`, and `MatShell`, with links to the official docs. If something on this page failed for you, that is a documentation bug worth reporting. Say which section, paste the command and the error, and attach the `--diagnose` output. --- ## Boundary interfaces Source: https://linstabpy.vercel.app/boundaries/ # Boundary interfaces This page is the reference for how you talk to LinStabPy from outside: CLI commands, YAML cases, Python imports, file inputs, and output files. Math conventions match the rest of the docs: $\omega = \mathrm{eig}(i L_0)$, unstable exactly when $\mathrm{Im}(\omega) > 0$, [resolvent](https://grokipedia.com/page/Resolvent_formalism) $H = C F (L_0 - (\gamma + i\omega) I)^{-1} F B$. ## Command line interface ### linstabpy (driver CLI) **Description**: Command-line entry for the shared case driver. It takes no case arguments. It reads the case path and output path from environment variables, checks the result cache by config hash and git SHA on rank 0, and either reuses cached results or runs the full driver (preflight, operator assembly, [PETSc](https://en.wikipedia.org/wiki/PETSc) and [SLEPc](https://en.wikipedia.org/wiki/SLEPc) solves, IO, provenance). **Source file**: `src/linstabpy/driver/cli.py` **Flags**: - `--version`: optional. Print the installed package version and exit without initializing PETSc or [MPI](https://grokipedia.com/page/Message_Passing_Interface). **Environment** (the real CLI surface, generated from `CLI_ENV_VARS`): - `LINSTABPY_CASE_CONFIG`: YAML case file. Aliases `LINSTABPY_CASE2_CONFIG`, `LINSTABPY_CASE1_CONFIG`. Default `./configs/parity_config.yaml`. - `LINSTABPY_CASE_OUTPUT_NPZ`: Result archive path. Aliases `LINSTABPY_CASE2_OUTPUT_NPZ`, `LINSTABPY_CASE1_OUTPUT_NPZ`. Default `{prefix}_{config_md5}_{git_sha8}.npz in the working directory`. - `LINSTABPY_CASE_LOG_PREFIX`: Log-line prefix and default output stem. Default `linstabpy`. - `LINSTABPY_FORCE_RECOMPUTE`: Set to 1, true, yes, or on to skip the cache check (CLI gate only; library runs always compute). Default `unset (cache check runs)`. - `LINSTABPY_STAGE_HEARTBEAT_SECONDS`: Progress heartbeat interval in float seconds. Default `30`. - `LINSTABPY_CASE_ROOT_ASSEMBLY`: Root-only assembly for iterative eigen ST solvers; stability-only. Aliases `LINSTABPY_CASE2_ROOT_ASSEMBLY`, `LINSTABPY_CASE1_ROOT_ASSEMBLY`. Default `unset (auto: on for iterative eigen ST solvers under MPI)`. !!! warning "Unknown arguments warn, then are ignored" The parser uses `parse_known_args`, so extra words on the command line do not error. Bare words and `--style` options print a `warning: ignoring unknown arguments [...]` line on stderr (rank 0 only) and the run continues with the env config; the exit code is unchanged. Single-dash PETSc options (`-ksp_monitor`, ...), plus one following value (`-ksp_rtol 1e-8`), still pass through to the solver stack. There is no `run` subcommand and no `--config` flag. An older revision of this page documented both. They never existed. Unknown `LINSTABPY_*` names in the environment get the same treatment: a stderr warning with a did-you-mean hint (registry: `src/linstabpy/config/env_registry.py`), then ignored. **Usage examples**: ```bash linstabpy --version ``` ```bash export LINSTABPY_CASE_CONFIG=path/to/case.yaml export LINSTABPY_CASE_OUTPUT_NPZ=results/out.npz mpiexec -n 8 linstabpy ``` ```bash python -m linstabpy.driver.cli ``` ### YAML case configuration interface **Description**: External configuration boundary for scientific cases. Users interact primarily by editing YAML case files that define baseflow paths, physics, numerics (including [finite-difference](https://grokipedia.com/page/Finite_difference_method) orders), boundary conditions, eigen/resolvent parameters, and solver policy. Files are ingested with a SafeLoader that accepts scientific floats, then structurally validated (TOP_LEVEL_ALLOWED keys such as paths/physics/numerics and related blocks). **Source File**: `src/linstabpy/utils/yaml_io.py; src/linstabpy/config/schema.py; src/linstabpy/utils/_solver_cfg.py` **Arguments**: - `config path` (str | Path): required - Filesystem path passed to load_config(path) returning a dict for driver/solvers. **Options**: - `paths`(mapping): optional - Top-level block for input/output paths (baseflow npy/npz dumps, mesh, artifacts). - `physics`(mapping): optional - Top-level physics parameters for the linearized [compressible](https://grokipedia.com/page/Compressible_flow) [Navier-Stokes](https://grokipedia.com/page/Navier–Stokes_equations) problem. - `numerics`(mapping): optional - Top-level numerics block including finite-difference order settings traversed by validation helpers. - `bc`(mapping): optional - Boundary recipe block: ordered `recipe` rules plus `penalty` (Dirichlet diagonal pin and Neumann row scale; driver default `1e8`, library and router default `1e4`). - `complex / int list fields`(complex | list[complex] | list[int]): optional - Loosely typed YAML scalars/lists parsed into Python complex and int lists (parse_complex, parse_complex_list, parse_int_list) for frequencies, [wavenumbers](https://grokipedia.com/page/Wavenumber), and similar solver inputs. **Top-level keys** (`TOP_LEVEL_ALLOWED`, generated — 73 keys): `L_cavity`, `Lx`, `Ly`, `Ma`, `Pr`, `Re`, `St`, `T_0`, `alpha_filter`, `alpha_temp`, `analysis`, `assembly_backend`, `baseflow_input`, `baseflow_preprocess`, `bc`, `bc_mask_layout`, `beta`, `cases`, `coarsen`, `conformal_scale`, `dedupe_tol`, `deform`, `dx`, `fd_order`, `floquet_shift`, `forcing_zero_vars`, `kappa`, `kx`, `m`, `mask_x_left`, `mask_x_right`, `model`, `n_modes`, `n_modes_resolvent`, `n_modes_stability`, `n_periods`, `ncv`, `numerics`, `nx`, `ny`, `omega`, `omega_shift_fallback`, `omega_target`, `paths`, `penalty`, `physics`, `preprocess`, `run_resolvent`, `run_spatial`, `run_stability`, `run_temporal`, `schema_version`, `solver`, `spatial_omega_imag_shift`, `sponge`, `stability_search_extra_modes`, `stride_eta`, `stride_ksi`, `temporal_quotas`, `temporal_targets`, `use_omega_shift_fallback`, `x_cavity_start`, `x_max_base`, `x_min_factor`, `x_periodicity`, `xrange`, `y_i`, `y_max`, `y_stretch_yi`, `y_symmetry`, `yrange`, `zero_dmudT`, `zero_dmudT_martini_parity` SBP `fd_order` values allowed: 2, 4, 6, 8. **Usage Examples**: ```bash load_config('cases/blasius_resolvent.yaml') ``` ```bash # case.yaml top-level keys constrained by validate_config / TOP_LEVEL_ALLOWED ``` ### make (tutorial docs build) **Description**: Secondary CLI boundary for building the tutorial PDF documentation with pdflatex (not part of the solver runtime). **Source File**: `docs/tutorial/Makefile` **Usage Examples**: ```bash make all ``` ```bash make clean ``` ```bash make tutorial.pdf ``` ## Integration Suggestions ### Python library import (lazy public API) Embed LinStabPy in larger openfluids or research pipelines by importing the package. Public symbols resolve on first attribute access so a plain import stays lightweight and does not eagerly initialize PETSc/MPI. Suitable for stability, [adjoint](https://grokipedia.com/page/Adjoint), [structural sensitivity](core-concepts.md), and resolvent workflows after the caller supplies baseflow data and config. **Example Code**: ``` import linstabpy # Heavy solvers load only on first attribute access print(linstabpy.__version__) # version without PETSc init # Configure via YAML dict or SolverConfig, then call driver/stability/resolvent APIs ``` **Public names** (`_LAZY_EXPORTS`, generated — 9 names): | Public name | Module | Attribute | |---|---|---| | `BaseFlow` | `linstabpy.baseflow` | `BaseFlow` | | `adjoint_eigenvalues` | `linstabpy.stability._eigenvalues` | `adjoint_eigenvalues` | | `build_lnse_operator` | `linstabpy.operators.assembly` | `build_lnse_operator` | | `build_lnse_system` | `linstabpy.operators.assembly` | `build_lnse_system` | | `compute_chu_norm` | `linstabpy.operators.energy_norm` | `compute_chu_norm` | | `resolvent_analysis` | `linstabpy.resolvent.analysis` | `resolvent_analysis` | | `resolvent_map` | `linstabpy.solvers.router.resolvent_map` | `resolvent_map` | | `temporal_eigenvalues` | `linstabpy.stability._eigenvalues` | `eigenvalues` | | `wavemaker` | `linstabpy.stability.sensitivity` | `wavemaker` | **Best Practices**: - Rely on lazy __getattr__ so import-time side effects stay minimal - Call PETSc/MPI initialization explicitly only when entering solver paths - Pass validated config dicts (validate_config) rather than ad-hoc nested structures - Prefer SolverConfig / from_config so YAML, env, and defaults share one precedence model ### YAML case file + driver orchestration Standard user-facing integration for batch and HPC runs: prepare mean-flow fields (npy/npz), write a Martini/Song-style YAML case, validate, and run through the driver CLI under MPI. Outputs are stability/resolvent artifacts, metrics, and provenance payloads. **Example Code**: ``` # case.yaml defines paths, physics, numerics, eigen/resolvent, solver policy from linstabpy.utils.yaml_io import load_config from linstabpy.config.schema import validate_config cfg = load_config('case.yaml') validate_config(cfg, source='case.yaml') # Then invoke driver.run / CLI entry with the merged config ``` **Best Practices**: - Keep dissipation apply flags aligned across stability and resolvent stages - Use schema validation before long MPI jobs to fail fast on unknown keys - Store baseflow and mesh paths relative to a reproducible case directory - Document expected artifacts alongside the case (module + NPZ name); the live example list is `run_all_examples.sh` `EXAMPLE_SPECS` ### PETSc / SLEPc / MUMPS / optional cuDSS backend External numerical stack integration. Core install uses NumPy/SciPy/PyYAML; extras enable HDF5, parallel PETSc/SLEPc (mpi4py, petsc4py, slepc4py), and optional GPU [cuDSS](https://docs.nvidia.com/cuda/cudss/). [MUMPS](https://en.wikipedia.org/wiki/MUMPS_(software)) ordering names are shared via mumps_options without importing PETSc at module load. **Example Code**: ``` # Install core + parallel stack # pip install 'linstabpy[petsc]' # or conda via environment.yml # Optional: pip install 'linstabpy[hdf5]' 'linstabpy[cudss]' # Runtime: export LINSTABPY_CASE_CONFIG=case.yaml && mpiexec -n N linstabpy ``` **Best Practices**: - Install matching mpi4py/petsc4py/slepc4py builds for the cluster MPI - Select MUMPS ordering (amd/pord/metis/qamd/auto) consistently in config - Treat cuDSS as optional GPU factor path shelled through PETSc PC - Do not import PETSc-dependent modules at package import time; follow lazy init patterns ### CFD baseflow data stores (file I/O) Ingest on-disk mean-flow fields as npy directories or npz archives (density, velocity components, temperature, mesh, sponge) produced by external nonlinear [CFD](https://grokipedia.com/page/Computational_fluid_dynamics). LinStabPy does not generate baseflows; integration starts at load/preprocess boundaries under baseflow/. **Example Code**: ``` # paths block in case.yaml points at baseflow archives # paths: # baseflow: ./data/meanflow.npz # output: ./artifacts/ ``` **Best Practices**: - Provide structured-mesh-compatible fields expected by grid_fields/preprocess - Apply sponge and shock-cleanup options via config rather than ad-hoc scripts - Version baseflow dumps alongside YAML for provenance --- ## Core concepts Source: https://linstabpy.vercel.app/core-concepts/ # Core concepts Short, exact definitions of the physics and math vocabulary used across this wiki. Each entry states what the thing is, the convention this repo uses, and where to read more. Solver objects (`KSP`, `PC`, `EPS`, `ST`, `MatShell`) live in the [solver glossary](solver-glossary.md). Code structure lives in the [module guide](modules/index.md). Shared conventions: the perturbation state is $q' = [\rho', u', v', w', T']^T$ per grid point; the assembled operator satisfies $\partial_t q' = L_0 q'$ in the interior; temporal eigenvalues are $\omega = \mathrm{eig}(i L_0)$ with $\exp(-i\omega t)$ time dependence. ## Stability vocabulary **Base flow.** A steady solution $\bar q(x, y)$ of the compressible [Navier-Stokes](https://grokipedia.com/page/Navier–Stokes_equations) equations, loaded from files and never computed here. Starting point of every analysis. Nondimensional, convective scaling: freestream density, velocity, and temperature equal 1. See [baseflow management](modules/baseflow-management.md). **Perturbation.** The small disturbance $q'(x, y, z, t)$ about the base flow. LinStabPy solves linear equations for $q'$ only. Nonlinear saturation is out of scope by construction. **Discrete temporal operator.** The assembled sparse matrix $iL_0$ whose [eigenvalues](https://grokipedia.com/page/Eigenvalues_and_eigenvectors) are the reported $\omega$. "Discrete" marks it as the matrix acting on grid vectors, as opposed to the continuous linearized operator. Every temporal eigensolve targets this matrix through shift-invert; the $i$ puts frequencies on the real axis and growth rates on the imaginary axis under the $\exp(-i\omega t)$ convention. **Temporal growth rate.** The imaginary part $\mathrm{Im}(\omega)$ of a temporal [eigenvalue](https://grokipedia.com/page/Eigenvalues_and_eigenvectors). Mode amplitude evolves as $\exp(\mathrm{Im}(\omega)\, t)$, so positive means exponential growth in time at a fixed station. Reported for every accepted mode with its postcheck residual. **Angular frequency.** The real part $\mathrm{Re}(\omega)$. Oscillations run as $\exp(-i\,\mathrm{Re}(\omega)\, t)$. [Shift-invert](solver-glossary.md) targets in YAML live in full complex $\omega$ space, so you aim at a frequency and growth rate together. **Direct mode.** A right eigenvector of the stability problem: the spatial shape $v(x, y)$ of the disturbance. Modes are normalized in the Chu norm (or $L^2$ where stated) before any comparison or plot. **Spatial amplification.** At prescribed real spanwise [wavenumber](https://grokipedia.com/page/Wavenumber) $\beta$ and real frequency $\omega$, the spatial problem returns a complex streamwise wavenumber $\alpha$. Amplification along $x$ is read from $\mathrm{Im}(\alpha)$ while following one branch; the sign interpretation depends on the branch, so never judge convective or absolute behavior from a single point. Build the spectrum first, then track the branch. See the tutorial and [stability analysis](modules/stability-analysis.md). **Generalized eigenvalue problem.** An eigenproblem of the form $A x = \lambda B x$ with $B \neq I$ (see [generalized eigenvectors](https://grokipedia.com/page/Generalized_eigenvector)). Ours is the spatial pencil $(L_c + \omega L_w)\, v = \alpha\, R\, v$. Boundary rows make $R$ singular, so $\lambda = \infty$ eigenvalues exist for the constrained subspace. Tiny or near-full problems go to dense [QZ](https://en.wikipedia.org/wiki/QZ_algorithm), which handles singular $B$ natively; production sizes use shift-invert [EPS](solver-glossary.md). **Boundary layer.** The thin region near a solid wall where viscosity brings the flow from freestream speed to rest at the wall. Blasius flow, the base flow of example 1, is the flat-plate similarity solution of the [boundary-layer](https://grokipedia.com/page/Boundary_layer) equations. Background: [Scholarpedia: Boundary Layer Theory](http://www.scholarpedia.org/article/Boundary_Layer_Theory). ## Resolvent vocabulary **Modal vs non-modal.** Modal stability asks about $t \to \infty$: do all eigenvalues sit on the stable side. Non-modal stability asks about finite time: how far can a perturbation wander before (possibly) decaying. For normal operators the answers agree; for non-normal ones they can differ wildly, and the finite-time answer is an optimization problem, not an eigenvalue. See Frantz, Loiseau and Robinet (2023), *Krylov methods for large-scale dynamical systems*, *Appl. Mech. Rev.* 75(3), 030802, modal and non-modal stability. **Resolvent.** The operator $H = C F (L_0 - (\gamma + i\omega) I)^{-1} F B$ mapping harmonic forcing to harmonic response at frequency $\omega$. Mathematically an instance of the [resolvent formalism](https://grokipedia.com/page/Resolvent_formalism). $B$ and $C$ restrict inputs and outputs (regions, variables), $F$ is an optional spatial filter inside the masks, $\gamma \ge 0$ is the discount below. Never formed as entries; applied through one factorization. The resolvent is the frequency-domain analogue of optimal-perturbation analysis below: same question (worst input, largest output), asked per frequency instead of per time horizon. See [resolvent analysis](modules/resolvent-analysis.md). **Optimal energy gain.** The largest achievable energy ratio $\sigma_1^2 = \max_{f \neq 0} \|Hf\|_E^2 / \|f\|_F^2$: response energy in the Chu norm over forcing energy in the forcing norm. Computed as $\sigma = \sqrt{\mathrm{Re}(\lambda)}$ from the weighted normal eigenproblem, never from an explicit [SVD](https://grokipedia.com/page/Singular_value_decomposition) of $H$. "Optimal" means optimal over all forcing shapes at that $(\beta, \omega)$, not over frequencies. **Forcing mode.** The optimal input shape $v$: the unit-norm vector attaining the leading gain. Normalized as $v^H W v = 1$ under the `chu` forcing norm or $v^H v = 1$ under `l2`. What to inject for maximum effect. **Response mode.** The optimal output shape $u = Hv/\sigma$: what grows when the optimal forcing is applied. Always recovered by applying $H$, which is why the forward residual sits at roundoff by construction. **Gain map.** Leading gains sampled over a $(\beta, \omega)$ grid: spanwise [wavenumber](https://grokipedia.com/page/Wavenumber) versus frequency. Peaks mark the flow's preferred amplifier settings. Production maps reuse anchor factorizations across nearby $\omega$ under a residual gate instead of refactoring per point. See [solver routing](modules/solver-routing.md). **Discount ($\gamma$).** A real shift added as $L_0 - (\gamma + i\omega) I$ with $\gamma > 0$ pushing poles off the imaginary axis. Keeps the shifted solve well posed near lightly damped eigenvalues. $\gamma = 0$ is the classical resolvent. Gains of a discounted resolvent answer the discounted question, so report $\gamma$ with the number. **Nonnormality.** $L_0$ does not commute with its [adjoint](https://grokipedia.com/page/Adjoint) in the energy inner product ($H^H W \neq W H$ in the weighted sense). Consequence: all modes stable does not mean quiet. Quantitatively, with $\lambda_r$ the leading real part and $\kappa(V) = \|V\|\|V^{-1}\|$ the eigenvector conditioning, $\exp(2\tau\lambda_r) \le \|\exp(\tau L)\|^2 \le \kappa(V)\,\exp(2\tau\lambda_r)$: the upper bound is where transient growth lives, and it grows with nonnormality. (In repo convention $\lambda_r = \mathrm{Im}(\omega)$.) Mechanically, picture almost-colinear eigenvectors interfering constructively before decaying; the extreme is a non-diagonalizable Jordan block with algebraic growth. See Frantz, Loiseau and Robinet (2023), *Krylov methods for large-scale dynamical systems*, *Appl. Mech. Rev.* 75(3), 030802, non-modal stability. Nonnormal operators support large transient amplification, which is exactly what the gain map measures. The related spectral picture is the [pseudospectrum](https://en.wikipedia.org/wiki/Pseudospectrum): $\{z : \|(zI - A)^{-1}\| > 1/\varepsilon\}$. Scholarpedia's [Stability](http://www.scholarpedia.org/article/Stability) article works the classic two-by-two transient-growth example. **Optimal perturbation and transient growth.** The finite-time gain $\mathcal{G}(\tau) = \max_{x_0} \|\exp(\tau L)x_0\|^2/\|x_0\|^2$ equals $\sigma_1^2$, the squared leading singular value of the propagator; the maximizer is the first right singular vector and the response is $\sigma_1$ times the first left one. With the adjoint available this recasts as an eigenproblem rather than an SVD. Textbook illustration: plane Couette flow at $Re = 300$, linearly stable yet with optimal gain near 100, streamwise vortices growing into streaks. See Frantz, Loiseau and Robinet (2023), *Krylov methods for large-scale dynamical systems*, *Appl. Mech. Rev.* 75(3), 030802, non-modal stability and illustration. **Receptivity.** The sensitivity of the flow to external disturbances: incoming noise, roughness, turbulence. Resolvent gains and [adjoint](https://grokipedia.com/page/Adjoint) modes quantify it. High gain at some $(\beta, \omega)$ means the environment can drive a large response there even when every eigenvalue is stable. Downstream of receptivity lies [transition to turbulence](http://www.scholarpedia.org/article/Transition_to_turbulence). ## Energy and adjoint vocabulary **Chu energy norm.** The compressible disturbance-energy inner product. Diagonal in our discretization: $W = \tfrac{1}{2}\,\mathrm{diag}(W_\rho, W_u, W_v, W_w, W_T)$ built from quadrature weights, base density and temperature, $\kappa$, and Mach number (times radius in axisymmetric mode). Defines "large" for gains, mode normalization, and adjoints. A zero diagonal entry is a hard error at shell construction, not a warning. **Adjoint mode.** An eigenvector of the energy-weighted adjoint $W^{-1} L_0^H W$. Where the direct mode shows what grows, the adjoint shows what the eigenvalue is sensitive to: the receptive shape. Returned in the temporal $\omega$ convention as conjugates of direct eigenvalues for real-valued base flows. **Structural sensitivity (wavemaker).** The pointwise product of direct and adjoint mode magnitudes (Giannetti-Luchini). Marks where a localized structural tweak moves the eigenvalue most: the wavemaker region. The control-oriented output of running direct plus adjoint analysis together. ## Numerics vocabulary **SBP operator.** A [summation-by-parts](https://grokipedia.com/page/Summation_by_parts) [finite-difference](https://grokipedia.com/page/Finite_difference_method) matrix: discrete derivative plus diagonal norm satisfying a discrete analogue of integration by parts, which bounds the discrete energy. Even interior orders 2, 4, 6, 8 only, each with validated boundary closures. Order 4 is the usual start. **Sponge layer.** An artificial damping region near outflow and farfield boundaries that absorbs outgoing waves before they reflect. Implemented as a damping field added to the operator during assembly. Related idea, different mechanism: the [perfectly matched layer](https://grokipedia.com/page/Perfectly_matched_layer). A zero sponge array means no damping, which is valid for some cases and wrong for others; the baseflow plots show it. **Boundary recipe.** The declarative list of [Dirichlet](https://grokipedia.com/page/Dirichlet_boundary_condition) and [Neumann](https://grokipedia.com/page/Neumann_boundary_condition) rules applied last in assembly: which borders, which variables, corner handling, curvilinear normals. Publishes the free degree-of-freedom sets every downstream solve must respect. **BC penalty.** The number behind `bc: penalty:`. For Dirichlet rows, assembly zeroes the row and column, then writes `penalty` on the diagonal, so the constrained equation reads `penalty` times $q_i = 0$: the DOF is pinned to zero and decoupled from interior equations. In the temporal EVP (on $iL_0$) these pinned DOFs sit near $i \times$ `penalty`, far from physical targets, and the free-DOF sets plus acceptance checks keep them out of results. For Neumann rows, the normal-derivative stencil row is multiplied by `penalty`: the homogeneous constraint is unchanged, but the row is scaled to balance against interior rows in the factorization. How to set it: keep the path default unless diagnostics complain. The driver and Song default to `1e8` (`bc.penalty` in YAML); assembly, the router, and the Martini examples default to `1e4` (nested `bc: penalty:` in their YAMLs, top-level `bc_penalty` or `penalty` in router dicts, `bc_penalty` kwarg in the direct API). Raise toward `1e6` (a validation case does exactly this) if constrained-Dirichlet leakage diagnostics show boundary modes leaking into results. Lower back down if you see conditioning trouble: a `1e8` diagonal next to order-1 interior entries is exactly the kind of scaling that can inflate residuals, and direct LU tolerates only so much. **Artificial dissipation.** The term $-\varepsilon Q_{\mathrm{diss}}$ subtracted from $L_0$, with $Q_{\mathrm{diss}}$ the Mattsson-Svard-Nordstrom operator. Damps grid-scale sawtooth modes the centered SBP stencils admit. Two halves must both be set (build the operator, set the strength) or nothing happens silently. Off on rectilinear meshes, on for meshes curved in both directions. See the Workflow page. **Implicit filter ($\alpha$).** The Mattsson-style spatial filter applied inside the resolvent masks and on the mesh. Strength parameter `alpha_filter`: $0.0$ is maximum filtering in the MATLAB convention, `"none"` disables. The adjoint path uses the adjoint-compliant variant so discrete adjointness holds. **Spanwise wavenumber ($\beta$).** The real [wavenumber](https://grokipedia.com/page/Wavenumber) of the $e^{i\beta z}$ spanwise dependence. Enters the linearized physics, so each $\beta$ needs its own assembled $L_0$. $\beta = 0$ is the two-dimensional problem. Gain maps sweep $\beta$ on the outer loop. **Floquet exponent.** The detuning parameter for spanwise-periodic [Floquet](https://grokipedia.com/page/Floquet_theory) analysis on streaky or modulated base flows. `continuum` is the default form; `matlab_parity` reproduces the LinStab2D algebra exactly for validation. Set per case, never guessed. **Matrix dual path.** Operators assemble locally in SciPy CSC (from CSR blocks) and convert to PETSc AIJ for distributed solves, or assemble AIJ directly where supported. Layouts, backends, and the conversion bridge are defined in the [solver glossary](solver-glossary.md). ## Workflow vocabulary **Shift target ($\sigma$).** The complex $\omega$-space point the eigensolver searches around. Lives in the YAML solver block, set explicitly on the spectral transformation. Near a branch it isolates modes cleanly; on a crowded real axis it returns mixtures. Moving slightly off-axis is a numerical search choice, not a physics change. Solver mechanics live in the [solver glossary](solver-glossary.md). **Acceptance.** The postchecks a computed mode must pass before it is trusted: full-operator residual, boundary-condition leakage, target distance, spectral boxes, conjugate completeness. Converged Ritz pairs fail these regularly on interior non-Hermitian problems, which is why acceptance exists as a separate stage from convergence. **Provenance.** The reproducibility record attached to every result: config hash, git SHA, backend and precision, mesh and norm settings. Stored inside the NPZ beside the numbers so a figure always traces back to the run that produced it. ## Beyond linear scope Linear analysis predicts onset. What happens after onset is nonlinear and outside this repo, but two landmarks orient the map. **Hopf bifurcation.** What happens after a mode goes unstable: an equilibrium sheds a limit cycle when a complex-conjugate eigenvalue pair crosses the imaginary axis. Linear analysis predicts the crossing ($\mathrm{Im}(\omega) = 0$ at criticality); the cycle itself needs the nonlinear terms. See [Scholarpedia: Andronov-Hopf bifurcation](http://www.scholarpedia.org/article/Andronov-Hopf_bifurcation) and [Bifurcation](http://www.scholarpedia.org/article/Bifurcation). **Transition to turbulence.** The full path from laminar flow to turbulence: receptivity, linear amplification (modal or transient), secondary instability, breakdown. This repo covers the linear amplification step and the sensitivity fields around it. See [Scholarpedia: Transition to turbulence](http://www.scholarpedia.org/article/Transition_to_turbulence). If a term appears elsewhere in this wiki without a definition or a link, that is a bug in that page. Link it here. --- ## FAQ Source: https://linstabpy.vercel.app/faq/ # Frequently asked questions Short answers to the questions new users actually ask, each pointing at the page with the full story. Conventions are fixed repo-wide: $\omega = \mathrm{eig}(i L_0)$, unstable exactly when $\mathrm{Im}(\omega) > 0$. ## Which analysis do I need? **Temporal stability** if you ask "does the flow sustain self-excited disturbances" — growth rates and frequencies of global modes. **Spatial stability** if you prescribe a real frequency and ask how a disturbance grows downstream — the unknown is the complex streamwise wavenumber $\alpha$ from a quadratic-turned-generalized problem of twice the size. **Adjoint modes** if you ask "where is the flow receptive" — sensitivity and wavemaker fields from $W^{-1}L_0^H W$. **Resolvent** if you ask "which external forcing is amplified most" — optimal gains with forcing and response shapes at real $\omega$. Temporal, adjoint, and resolvent act on the same assembled $L_0$; spatial builds its own pencil from the same coefficients. Details in [core concepts](core-concepts.md), solver mechanics in the [solver glossary](solver-glossary.md). ## My spectrum moved when I changed the shift. Is that normal? Yes, if you changed which branch the solver locked onto. Shift-invert finds eigenvalues near the target; on a crowded real axis, neighboring targets can return mixtures from different branches. Fix the branch first (small `nev`, a few nearby targets, compare postcheck residuals and target distances), then scale up. Moving the target slightly off the real axis often isolates a branch more cleanly. This is a numerical search choice, not a physics change. See [getting started](getting-started/index.md#your-first-own-case) and the paper's note that slow convergence near branches means widening the subspace or moving the shift. ## Do I turn the dissipation on? Run the grid-scale leak test first: feed the operator a field alternating along one grid direction and constant along the other; a transverse derivative must return zero, so $\|D_x q\|/\|D_y q\|$ measures pure metric mixing. At round-off (the Song and Savarino meshes read below $10^{-11}$), the mesh is rectilinear: leave $\varepsilon = 0$ and confirm insensitivity under refinement. At order one (0.26 on one LES curved-wall mesh, 0.10–0.12 on a wavy-wall test), the mesh is bent: switch dissipation on and tune $\varepsilon$ on that mesh until the gain or eigenvalue stops moving. A value tuned elsewhere does not transfer. And remember the two halves: building the operator and setting the strength are separate flags, and setting only one is a silent no-op. Costs and the full procedure come from the paper's sawtooth section: $\varepsilon = 10^{-2}$ damps grid content at rate $\approx 0.3$ while resolved scales see $\varepsilon \times 10^{-2}$, for 15% more nonzeros and 3–4× wall time. ## My resolvent response shows point-to-point wiggles. What is that? The odd–even (sawtooth) mode: the $q_i = (-1)^i$ family that every centered first-derivative stencil annihilates exactly, at any order. Only viscous terms damp it, at rate $\sim 1/(\mathrm{Re}\,h^2)$, so at high $\mathrm{Re}\,h^2$ it sits nearly neutral and whatever projects onto it stays. On rectilinear meshes it decouples and refinement removes it; on bent meshes the cross-metrics couple it into the resolved field and you need the dissipation from the previous answer. See the sawtooth discussion in [core concepts](core-concepts.md#numerics-vocabulary). ## What does `bc: penalty:` do, and what value? Dirichlet rows are zeroed and the diagonal is set to `penalty`, pinning the DOF (`penalty`·$q_i = 0$); pinned modes land near $i \times$ `penalty`, far from physical targets. Neumann rows are scaled by it for row balance. Driver default `1e8`, library and router default `1e4`. Keep the default unless constrained-leakage diagnostics complain, then raise toward `1e6` first. Full mechanics in [core concepts](core-concepts.md#numerics-vocabulary). ## Why does `linstabpy run case.yaml` run the wrong case? Because that command never existed. The CLI takes no case arguments — only `--version`. Case and output travel by `LINSTABPY_CASE_CONFIG` and `LINSTABPY_CASE_OUTPUT_NPZ`, and unknown bare words print a stderr warning but are otherwise ignored by design of the parser. See [getting started](getting-started/index.md#the-cli-as-it-really-is). ## Why did my run print `Cache valid` and finish instantly? Rank 0 found an output NPZ whose stored config hash and git SHA match the current ones, so it skipped the compute. Point each run at its own output path, or export `LINSTABPY_FORCE_RECOMPUTE=1` deliberately. See [getting started](getting-started/index.md#the-cli-as-it-really-is). ## Which SBP order, and why only even ones? Order 4 to start. Only even orders 2, 4, 6, 8 ship because only they have validated diagonal-norm boundary closures in this family. Higher order buys accuracy per point but stiffens the operator and grows stencils. See [core concepts](core-concepts.md#numerics-vocabulary). ## How big should `ncv` be? Bigger than `nev`, then bigger again near branches. Measured guidance from the paper's campaign: on a 450k-unknown resolvent point, cutting the Krylov space from 20 vectors to 5 cut solves per point from 40 to 10 and moved the gain by less than $10^{-10}$ — subspace size sets cost, not precision, once past the threshold. If accepted-mode quality plateaus while the solve stays active, the subspace is too small for the neighborhood, not the tolerance too loose. Defaults per path live in the [solver glossary](solver-glossary.md#slepc-eps-eigensolvers-and-st-spectral-transformations). ## `nconv > 0`. Am I done? No. Converged Ritz pairs can still fail full-operator residuals, boundary leakage, or target distance on interior non-Hermitian problems. Trust the acceptance postchecks, not the convergence count: requested `nev` vs converged `nconv`, explicit postcheck residuals, leakage diagnostics, and target distances for accepted and rejected modes alike. See [getting started](getting-started/index.md#your-first-own-case). ## Can I compare my gain with another paper's number? Only if both weights match: response weight $W$ and forcing weight $W_f$. Kinetic-energy-only norms misrank compressible modes because thermodynamic fluctuations carry real energy share, and a Chu gain is simply not the same number as an $L^2$ gain. Always report forcing norm, response norm, $\beta$, $\omega$, $\gamma$, grid, and SBP order with the value. See [core concepts](core-concepts.md#resolvent-vocabulary). ## Direct or iterative linear solves? Direct (`KSPPREONLY` + `PCLU` + MUMPS) unless memory says otherwise — it is the accuracy reference, and one factorization serves every iteration and every $H$/$H^H$ apply at that shift. Iterative paths (`fgmres_lu`, `strumpack_fgmres`) are explicit opt-ins for cases that do not fit direct factorization, always with `true_residual` on because transformed residuals can lie. The decision procedure with measured costs is in the [solver glossary](solver-glossary.md#choosing-a-configuration). ## One rank or many? Factorization dominates (about four fifths of a resolvent point: 60 s factor vs 0.44 s solves and 5.8 s eigensolve at 450k unknowns), and it scales poorly: measured 190.5 s at 4 ranks, 154.6 s at 8, 148.2 s at 16, while peak memory doubled from 18.6 to 40.6 GiB. Start at 4 ranks with one thread each, compare on your machine, and scale out across cases (scheduler arrays) before scaling up ranks. Numbers from `examples/savarino/*.csv`; protocol in [getting started](getting-started/index.md#the-example-family-map). ## Where do I put a question this page does not answer? Open an issue with the section that failed you, the command, the error, and the `--diagnose` output. Install problems especially — they are the most useful reports this project gets. --- ## Overview Source: https://linstabpy.vercel.app/overview/ # Project overview | Attribute | Value | |-----------|--------| | **System** | LinStabPy | | **Project type** | CLI and scientific library | | **Domain** | Compressible flow linear stability and resolvent analysis | | **Audience** | Master students and researchers fluent in math, linear algebra, and HPC, plus CFD engineers and scientific developers | --- !!! info "Conventions" The discrete perturbation state is $q' = [\rho', u', v', w', T']^T$ at each grid point. Assembly builds the [sparse](https://grokipedia.com/page/Sparse_matrix) matrix $L_0$ with $\partial_t q' = L_0 q'$ in the interior, then applies sponge, dissipation, and boundary conditions in that fixed order. Temporal [eigenvalues](https://grokipedia.com/page/Eigenvalues_and_eigenvectors) are $\omega = \mathrm{eig}(i L_0)$ with time factor $\exp(-i\omega t)$. Growth in time is $\mathrm{Im}(\omega) > 0$. Targets you pass in YAML live in $\omega$ space. This matches `stability/_eigenvalues.py`, which is the source of truth. The [resolvent](https://grokipedia.com/page/Resolvent_formalism) at $(\beta, \omega)$ with real discount $\gamma$ is $H = C F (L_0 - (\gamma + i\omega) I)^{-1} F B$. Gains are weighted [singular values](https://grokipedia.com/page/Singular_value_decomposition). The default energy is the diagonal Chu norm described in Architecture. [SBP](https://grokipedia.com/page/Summation_by_parts) orders are exactly $2$, $4$, $6$, $8$ (even only, with validated boundary closures). The default shifted solve is `KSPPREONLY` with `PCLU` through [MUMPS](https://en.wikipedia.org/wiki/MUMPS_(software)), driven by [SLEPc](https://en.wikipedia.org/wiki/SLEPc) `EPS` with Krylov-Schur and shift-invert. Every solver term here is defined in the [solver glossary](solver-glossary.md). ## Project introduction ### Project name and description **LinStabPy** is a Python toolkit for **linear stability** and **resolvent analysis** of [compressible](https://grokipedia.com/page/Compressible_flow) flows. It started as a Python port of LinStab2D style MATLAB workflows and now runs distributed on HPC through [MPI](https://grokipedia.com/page/Message_Passing_Interface), [PETSc](https://en.wikipedia.org/wiki/PETSc), and SLEPc. A run starts from two inputs you supply. One is the **base flow**, a steady compressible [Navier-Stokes](https://grokipedia.com/page/Navier–Stokes_equations) solution stored as arrays on a structured grid. The code linearizes about it and never computes it. The other is a short **YAML case file** with physics, numerics, and solver policy. From there the pipeline builds the **linearized Navier-Stokes (LNSE)** matrix $L_0$ with **summation-by-parts (SBP)** [finite differences](https://grokipedia.com/page/Finite_difference_method), enforces boundary conditions, and solves one or more of: * **Temporal** eigenproblem, $\omega = \mathrm{eig}(iL_0)$ near a target * **Spatial** eigenproblem, generalized $Ax = \lambda Rx$ for streamwise [wavenumber](https://grokipedia.com/page/Wavenumber) * **[Adjoint](https://grokipedia.com/page/Adjoint)** eigenproblem, $W^{-1} L_0^H W$ for receptivity and sensitivity * **Resolvent** problem, leading [gains](core-concepts.md) and forcing and [response modes](core-concepts.md) of $H$ Large sparse algebra runs through **PETSc and SLEPc**. The spine is **factorize once, solve many**. One sparse direct factorization of the shifted operator serves all [Krylov](https://grokipedia.com/page/Krylov_subspace) iterations and all $H$ and $H^H$ applies at that shift. The default backend is **MUMPS**. An optional **[cuDSS](https://docs.nvidia.com/cuda/cudss/)** path exists for GPU factors. ### What you get | Capability | What it gives you | |------------|-------------------| | Large LNSE eigen and resolvent studies on HPC without MATLAB | Same analysis semantics as LinStab2D, on clusters, with a reproducible Python stack | | Mean flow files turned into modes, gains, and sensitivity fields | A direct bridge from RANS, LES, or DNS averages to transition and control insight | | Temporal, spatial, and adjoint spectra with structural sensitivity | Growth rates, frequencies, wavenumbers, and wavemaker maps | | Resolvent gain maps over frequency $\omega$ and spanwise wavenumber $\beta$ | Quantified nonnormal amplification, even where all modes are stable | | YAML driven runs with schema checks and provenance | Config hash, git SHA, backend, and cost recorded with every result | In short, LinStabPy turns [CFD](https://grokipedia.com/page/Computational_fluid_dynamics) mean flow dumps into stability spectra, resolvent gains, and sensitivity results for **[transition](http://www.scholarpedia.org/article/Transition_to_turbulence)** and **flow control** work. ### How it is built * **Language and delivery.** Python CLI plus importable library. Plain `import linstabpy` stays light and does not initialize PETSc or MPI. The driver initializes them when a solve stage starts. * **Discretization.** Structured grids. SBP differences at even orders $2$ to $8$. Cartesian and cylindrical (axisymmetric) coefficient paths. Optional Mattsson style filters and artificial dissipation. * **State and operators.** LNSE for $q' = [\rho', u', v', w', T']$. Sparse $L_0$ with sponge and [Dirichlet](https://grokipedia.com/page/Dirichlet_boundary_condition) and [Neumann](https://grokipedia.com/page/Neumann_boundary_condition) recipes applied last. * **Solvers.** SLEPc `EPS` for eigenproblems. Resolvent through a weighted normal `MatShell` with repeated shifted solves. * **Parallelism.** MPI distribution of matrices, vectors, and case orchestration through `petsc4py` and `slepc4py`. Complex scalars are required for production paths. * **Linear algebra core.** PETSc `KSP` and `PC`. Sparse direct [LU](https://grokipedia.com/page/LU_decomposition) through MUMPS by default. Optional GPU sparse direct through cuDSS as a PETSc [preconditioner](https://grokipedia.com/page/Preconditioner) shell. * **Local numerics and IO.** NumPy and SciPy for arrays and CSC and CSR assembly before PETSc conversion. Base flow from `npy` directories or `npz` archives. Cases from YAML. * **Architecture style.** Staged pipeline with a dedicated **solver router** for cost aware $\beta$ and frequency campaigns with anchors, factor reuse, and residual gates. ## Who uses it ### Roles ```mermaid flowchart LR subgraph users [Users] FDR["Fluid dynamics researchers"] ENG["CFD and aerospace engineers"] DEV["Scientific software developers"] end SYS["LinStabPy"] FDR -->|"spectra, adjoints, sensitivity, maps"| SYS ENG -->|"base flow ingest, BC recipes, sweeps"| SYS DEV -->|"library API, typed config, provenance"| SYS ``` | Role | Description | |------|-------------| | :microscope: **Fluid dynamics researchers** | People studying hydrodynamic stability, transition, and receptivity in compressible flow. | | :small_airplane: **CFD and aerospace engineers** | People who already hold a mean flow and need its linear modes and amplifier behavior. | | :material-laptop: **Scientific software developers** | People who embed LNSE stages in larger pipelines and need a clean API and contracts. | If you are a master student, start in the researcher row. You need the base flow idea, the $\omega$ convention above, and one example. MPI and Krylov details can wait until your first map campaign feels slow. ### Usage scenarios | Scenario | Who | What happens | |----------|-----|--------------| | **Temporal eigenvalue analysis** | Researchers, engineers | Load mean flow and mesh, assemble $L_0$ with boundary conditions, factorize the shifted system, run SLEPc near the target, emit $\omega$, modes, residuals, and provenance. | | **Resolvent gain analysis** | Researchers, engineers | Build $H$ from $L_0$ at target $\omega$, solve the weighted normal eigenproblem, return gains with forcing and response pairs. | | **Gain map over $\beta$ and $\omega$** | Researchers on HPC | One $L_0$ per spanwise wavenumber $\beta$, frequency bands marked, anchors factorized, nearby $\omega$ points reuse factors as preconditioners under a residual gate. | | **Spatial eigenproblem** | Researchers, engineers | Assemble spatial block pair $(A, R)$, run generalized SLEPc for streamwise wavenumber and spatial growth. | | **Base flow ingest and operator prep** | Engineers, developers | `npy` or `npz` dumps pass through sponge, shock cleanup, and viscosity setup, then mesh and SBP, then BC ready $L_0$. | | **Library embedding** | Developers | Import analysis functions without PETSc init on import, use typed config, SciPy and PETSc bridges, metrics and provenance. | ### What each role needs **Researchers need:** * Temporal and spatial spectra with a clear $\omega$ convention * Adjoint modes and [structural sensitivity](core-concepts.md) ([wavemaker](core-concepts.md)) fields * [Gain maps](core-concepts.md) over $\omega$ and $\beta$ with residuals attached * YAML cases that rerun bit reproducibly on MPI clusters **Engineers need:** * Ingest of `npy` and `npz` base flows without reformatting by hand * Dirichlet and Neumann recipes expressed as data, not code edits * Direct sparse solves with MUMPS, with GPU factors available * Sweeps that reuse factors instead of refactoring at every point **Developers need:** * Importable API that stays free of PETSc init on import * Typed config with schema validation * SciPy and PETSc conversion helpers in one place * Provenance and metrics with stable key names ## System boundaries ### Scope **In scope** starts from **loaded mean flow and mesh geometry** and ends at **results, metrics, and payloads**. Inside that envelope, LinStabPy: 1. Assembles $L_0$ and spatial pairs on structured SBP grids 2. Enforces boundary condition recipes 3. Factorizes shifted systems once per shift 4. Runs eigen and resolvent workflows, including parametric maps 5. Orchestrates cases through CLI and driver code with validation, IO, and provenance It does **not** generate nonlinear base flows or ship visualization products. ```mermaid flowchart TB subgraph outside [Outside the boundary] CFD["Nonlinear CFD solvers"] GUI["Plotting and GUIs"] UNS["Unstructured or full 3D frameworks"] MULTI["General non-LNSE multiphysics"] CLOUD["Cloud auth and tenancy"] end subgraph inside [LinStabPy] BF["Base flow load and preprocess"] MS["Structured mesh, SBP, filters"] OP["LNSE coefficients and assembly"] BC["BC recipes"] ST["Temporal, spatial, adjoint eigensolvers"] RS["Resolvent gains and maps"] FA["Sparse factorization"] RT["Solver routing"] DR["CLI, config, IO, provenance"] end CFD -.->|"mean flow files only"| BF BF --> MS --> OP --> BC OP --> ST OP --> RS FA --> ST FA --> RS RT --> RS DR --> BF DR --> ST DR --> RS DR --> RT ``` ### Included components | Component | Role | |-----------|------| | Base flow load and preprocess | Ingest mean fields and make them safe to linearize | | Structured mesh with SBP and filters | Discrete derivatives, norms, quadrature, stabilization | | LNSE coefficients and assembly | Build sparse $L_0$ and spatial pairs | | Boundary condition recipes | Enforce Dirichlet and Neumann rows and publish free degrees of freedom | | Temporal, spatial, and adjoint eigensolvers | Core stability product | | Resolvent operator, gains, and maps | Amplifier and receptivity product | | Sparse factorization layer | Factorize once, solve many | | Solver routing with anchors and sweeps | Cost aware parametric campaigns | | CLI and driver with config, IO, and provenance | Reproducible HPC entry and exit | **Internal modules.** Case Orchestration, Configuration Management, Baseflow Management, Mesh Discretization, LNSE Operator Assembly (including boundary conditions), Stability Analysis, Resolvent Analysis, Solver Routing and Sweep Optimization, Linear Algebra Infrastructure, Shared Scientific Utilities. ### Excluded components | Excluded | Why | |----------|-----| | Nonlinear CFD solvers | Upstream data providers only, coupled through files | | GUIs and interactive visualization | Results are files for external tools | | Unstructured or fully 3D volume frameworks | Scope is structured grids | | General non-LNSE multiphysics | Single purpose toolkit | | Cloud deployment, auth, and tenancy | HPC CLI and library, not a hosted service | ## External pieces ### List | External piece | Type | What it does here | |----------------|------|-------------------| | **PETSc** | Library | Distributed `Mat` and `Vec`, `KSP` and `PC`, parallel runtime | | **SLEPc** | Library | `EPS` eigensolvers for temporal, spatial, adjoint, and resolvent normal problems | | **MUMPS** | Backend via PETSc | Default sparse direct LU for shifted systems | | **MPI** | Runtime | Process parallelism for matrices, vectors, and orchestration | | **cuDSS** | Optional backend | Optional GPU sparse direct path, wrapped as a PETSc shell | | **NumPy and SciPy** | Libraries | Arrays, CSC and CSR assembly, local kernels before PETSc conversion | | **Base flow file stores** | File input | Mean fields as `npy` directories or `npz` archives | | **YAML case config** | File input | Case, mesh, eigen, resolvent, and solver policy | ### How they connect ```mermaid flowchart LR subgraph inputs [File inputs] YAML["YAML case config"] CFDData["Base flow files"] end subgraph linstab [LinStabPy] Core["LNSE stability and resolvent"] end subgraph libs [Libraries and runtimes] NP["NumPy and SciPy"] MPI["MPI"] PETSc["PETSc"] SLEPc["SLEPc"] MUMPS["MUMPS"] cuDSS["cuDSS, optional"] end YAML -->|"merge and validate"| Core CFDData -->|"mean flow and mesh"| Core Core -->|"local arrays, CSC and CSR"| NP Core -->|"distributed runtime"| MPI Core -->|"Mats, Vecs, KSP, PC"| PETSc PETSc --> MUMPS PETSc -.-> cuDSS Core -->|"EPS problems"| SLEPc SLEPc --> PETSc ``` **File inputs.** * **YAML case config.** Merged with environment overrides and defaults, then schema checked for eigen, resolvent, mesh, MPI, and router policy before any heavy stage runs. * **Base flow file stores.** Nondimensional primitives and mesh plus optional sponge fields. No live coupling to a CFD process. **Libraries and runtimes.** * **NumPy and SciPy.** Host side fields and sparse assembly before distributed conversion. * **MPI with PETSc.** Distributed execution fabric, matrices and vectors, direct and preconditioned solves. * **MUMPS through PETSc.** Primary sparse direct backend. * **cuDSS, optional.** GPU direct factor where selected by resource policy. * **SLEPc.** Eigenvalue and gain eigenproblems on PETSc objects. ### Which dependencies matter most | Dependency | Criticality | Notes | |------------|-------------|-------| | PETSc with MPI | Required for distributed paths | Core parallel algebra and orchestration runtime | | SLEPc | Required for eigen and resolvent outputs | Built on PETSc structures | | MUMPS | Default direct path | Reached as a PETSc backend, not a separate app | | NumPy and SciPy | Required | Local numerics and IO bridge | | cuDSS | Optional | Probed for capability, chosen under RAM and cost policy | | Base flow files | Required as data source | One way file dependency | | YAML configs | Required for CLI campaigns | Policy source of truth after merge and validation | LinStabPy **consumes** CFD artifacts and HPC libraries. Edges are library linkage and file IO, not service APIs. ## System context ### Context diagram ```mermaid graph TB Researcher["Fluid dynamics researcher"] Engineer["CFD and aerospace engineer"] Developer["Scientific software developer"] LinStabPy["LinStabPy: compressible LNSE stability and resolvent on HPC"] YAMLConfig["YAML case config"] CFDStore["Base flow file stores"] NumPySciPy["NumPy and SciPy"] MPIRuntime["MPI"] PETScLib["PETSc"] SLEPcLib["SLEPc"] MUMPSLib["MUMPS"] CuDSSLib["cuDSS, optional"] Researcher -->|"YAML cases, MPI jobs, spectra, gains"| LinStabPy Engineer -->|"base flows, BC policy, sweeps"| LinStabPy Developer -->|"API, typed config, metrics"| LinStabPy YAMLConfig -->|"case and solver policy"| LinStabPy CFDStore -->|"mean flow and mesh"| LinStabPy LinStabPy -->|"arrays, CSC and CSR"| NumPySciPy LinStabPy -->|"process parallelism"| MPIRuntime LinStabPy -->|"distributed Mats, Vecs, KSP, PC"| PETScLib LinStabPy -->|"eigen and gain problems"| SLEPcLib PETScLib -->|"sparse direct LU"| MUMPSLib PETScLib -.->|"optional GPU factor"| CuDSSLib SLEPcLib -->|"uses PETSc objects"| PETScLib ``` ### Main flows #### Flow A. Temporal eigenvalue analysis 1. Launch CLI or call the stability API. 2. Initialize PETSc and MPI. Merge and validate YAML, env, and defaults. 3. Load and preprocess mean flow. Build structured mesh with SBP. 4. Compute LNSE coefficients. Assemble sparse $L_0$. Apply BC recipes. 5. Factorize the shifted operator at the $\omega$ target. 6. Run SLEPc near the target. Validate modes by residuals and boundary leakage. 7. Emit eigenvalues, vectors, metrics, and provenance. #### Flow B. Resolvent gain analysis 1. Resolve frequency, weights, $B$ and $C$ maps, and backend policy. 2. Reuse the same base flow, mesh, and $L_0$ foundation as Flow A. 3. Factorize $(L_0 - (\gamma + i\omega) I)$. Apply $H$ and $H^H$ through that factor. 4. Run SLEPc on the weighted normal shell. Extract gains with $\sigma = \sqrt{\mathrm{Re}(\lambda)}$. 5. Package gains, modes, residuals, and solver metrics. #### Flow C. Gain map over $\beta$ and $\omega$ 1. Validate map, router, and resource policy. 2. Assemble one $L_0$ per spanwise wavenumber $\beta$. 3. Mark frequency bands. Walk lines with **anchor factorizations**, preconditioner reuse, residual gated refactor, and optional stable sweeps. 4. Probe MUMPS, LU, and cuDSS backends under RAM and cost policy. 5. Record leading gains with cost samples and certification metadata. #### Flow D. Spatial stability 1. Resolve `solver.spatial` policy, with fallback to `solver.eigen`. 2. Assemble spatial block pair from prepared coefficients and mesh. 3. Configure direct shift-invert backend. Run generalized SLEPc solve for streamwise wavenumber and growth. ### Decisions visible at this level | Decision | What it means | Consequence | |----------|---------------|-------------| | **Port MATLAB workflows to Python HPC** | Keep LinStab2D analysis semantics, add MPI orchestration | Familiar quantities, no MATLAB runtime | | **File based CFD boundary** | Mean flow enters only as `npy` and `npz` | Clean split from nonlinear solvers | | **PETSc and SLEPc as algebra stack** | No custom distributed eigensolver | Strong fit with cluster numerical ecosystems | | **Factorize once, solve many** | Direct factors serve shift-invert and resolvent applies | MUMPS and optional cuDSS appear explicitly | | **YAML with schema as contract** | Validate before heavy compute | Config files are first class inputs | | **No GUI and no cloud tenancy** | CLI with library and file payloads | Plotting and multi tenant ops stay outside | | **Structured SBP only** | Exclude unstructured and full 3D frameworks | Narrow, deep scope | | **Lazy library import** | Public API without PETSc init on import | Embeddable without forcing MPI at import time | ## Technical architecture ### Stack | Layer | Choice | |-------|----------------| | **Language and packaging** | Python, CLI entry plus importable package | | **Configuration** | YAML cases, env overrides, defaults, typed models, schema checks | | **Local arrays and sparse** | NumPy, SciPy CSC and CSR | | **Discretization** | Structured mesh, SBP orders $2$, $4$, $6$, $8$, filters, pole treatments | | **Parallel runtime** | MPI | | **Distributed algebra** | PETSc `Mat` and `Vec`, `KSP`, `PC` | | **Eigen and gain algebra** | SLEPc `EPS` | | **Sparse direct backends** | MUMPS by default through PETSc, optional cuDSS | | **Data in** | `npy` directories, `npz` archives, YAML | | **Data out** | Result payloads with eigenvalues and vectors, gains, modes, metrics, provenance | ### Patterns | Pattern | Where it shows up | |---------|-------------------| | **Staged pipeline** | Config, base flow, mesh, operators and BCs, factorize, solve, payload | | **File and library boundary** | Domain science inside, CFD files and PETSc and SLEPc and MPI outside | | **Factorize once, solve many** | Shared layer under shift-invert eigen and resolvent applies | | **Strategy and backend registry** | MUMPS, LU, and cuDSS selection under cost and RAM policy | | **Router optimized sweep** | Anchors, recycling, residual gates, stable sweeps for $\beta$ and $\omega$ maps | | **Declarative BC recipes** | Dirichlet and Neumann rules sequenced onto operators with free degree sets | | **Lazy facade API** | Public surface defers heavy runtime init | ### Key design choices 1. **Single policy source.** Merged, validated YAML with env and defaults drives mesh, eigen, resolvent, MPI, and router behavior. 2. **Single operator family.** Stability, adjoint, sensitivity, and resolvent consume BC enforced assembled operators, not ad hoc matrices per product. 3. **HPC first.** MPI distributed vectors and matrices with case orchestration are first class. MATLAB parity is a workflow goal, not the scaling ceiling. 4. **Cost aware maps.** Large $\beta$ and $\omega$ studies use the router. Factorization cost is explicit in metrics and certification, not hidden in a double loop. 5. **Strict boundary.** Nonlinear CFD, visualization, unstructured 3D, general multiphysics, and cloud tenancy stay out of scope. 6. **Dual use.** CLI and driver for full cases, importable library for embedding. Both share the same pipeline without forcing PETSc on import. ### What connects most strongly inside At context level the software is one bubble. The strongest internal relations, which explain why PETSc and SLEPc and MUMPS matter, are: * Case Orchestration calls Configuration, Baseflow, Stability, Resolvent, Routing, Operators, and Linear Algebra * Operator Assembly composes Mesh and consumes Baseflow data * Stability and Resolvent consume Operators for $L_0$ and Linear Algebra for factorized shifted solves * Routing calls Resolvent and Linear Algebra for map walks and anchor reuse Package detail lives in Architecture and the Module guide. ## Summary **LinStabPy** is an HPC Python **CLI and library** for **compressible LNSE stability and resolvent analysis** on **structured SBP grids**. It starts from **base flow files** and **YAML cases** and delivers **spectra, modes, gains, sensitivity, metrics, and provenance**. It excludes CFD generation, GUIs, and cloud service layers. **People** are researchers, engineers, and developers. **System** is LinStabPy. **External pieces** are YAML and base flow file stores, NumPy and SciPy, MPI, PETSc, SLEPc, MUMPS, and optional cuDSS. This overview sets the boundary for the rest of the docs. --- ## Architecture Source: https://linstabpy.vercel.app/architecture/ # Architecture **LinStabPy** is a Python CLI and library for linear stability and resolvent analysis of compressible viscous flows. It ports LinStab2D style linearized Navier-Stokes (LNSE) workflows from MATLAB to an MPI enabled HPC stack on [PETSc](https://en.wikipedia.org/wiki/PETSc) and [SLEPc](https://en.wikipedia.org/wiki/SLEPc), with sparse direct factorization ([MUMPS](https://en.wikipedia.org/wiki/MUMPS_(software)) by default, optional [cuDSS](https://docs.nvidia.com/cuda/cudss/)). This page describes the packages, how they connect, the main processes, and how a case runs. It follows the tree and the module notes in `docs/architecture.md`. !!! info "Conventions" Temporal [eigenvalues](https://grokipedia.com/page/Eigenvalues_and_eigenvectors) are $\omega = \mathrm{eig}(i L_0)$ with $\exp(-i\omega t)$ time dependence. Unstable means $\mathrm{Im}(\omega) > 0$. The shift-invert target lives in $\omega$ space. The [resolvent](https://grokipedia.com/page/Resolvent_formalism) is $H = C F (L_0 - (\gamma + i\omega) I)^{-1} F B$ with $\sigma = \sqrt{\mathrm{Re}(\lambda)}$ from the weighted normal problem. [SBP](https://grokipedia.com/page/Summation_by_parts) orders are $2$, $4$, $6$, $8$ only. The default shifted solve is `KSPPREONLY` with `PCLU` through MUMPS, driven by SLEPc Krylov-Schur with `ST=SINVERT`. Assembly order is fixed: raw blocks, `inv(B/1j)` scaling, sponge, dissipation, then BC rows last. ## Architecture overview ### Architecture design philosophy LinStabPy is designed as a **staged scientific pipeline**, not a service mesh or multi-tenant application. The system starts from [CFD](https://grokipedia.com/page/Computational_fluid_dynamics) mean-flow dumps and YAML case policy, and ends at spectra, modes, [gains](core-concepts.md), metrics, and provenance payloads. Core design tenets: | Principle | Manifestation in the codebase | |-----------|-------------------------------| | **Staged pipeline** | Driver coordinates discrete stages: config → baseflow → mesh → L0 → factorize → solve → payload, with stage timing and heartbeats | | **Factorize-once, solve-many** | Shared `linalg.factorize` core underpins shift-invert eigenproblems and resolvent \(H\) applies | | **Config as single policy source** | YAML + environment + defaults → schema validation → typed `Config` subsections | | **Lazy loading** | PEP 562 `__getattr__` on package, mesh, BC, linalg, stability, resolvent, and router so help/version paths stay light and do not initialize PETSc/MPI | | **SciPy and PETSc dual paths** | Local CSC/CSR assembly is possible; production path converts to distributed PETSc `Mat` | | **Separation of concerns** | Physics (coefficients/assembly) vs. numerics (mesh/SBP) vs. routing (cost-aware sweeps) | | **Reproducibility** | Config hash, git SHA, metrics schema, fixed NPZ contract keys | The architectural style is a **modular layered scientific toolkit**: ```text Orchestration → Domain Analysis → Discretization / Operators → Linear Algebra → Utilities ``` ### Core architecture patterns 1. **Pipeline and staged orchestration.** `driver/run.py` sequences the scientific stages in a fixed order. 2. **Shared kernel with multiple products.** One BC enforced LNSE operator family feeds temporal eigen, [adjoint](https://grokipedia.com/page/Adjoint) and sensitivity, spatial eigen, single frequency resolvent, and $\beta$ and $\omega$ maps. 3. **Strategy with backend registry.** Factor backends (MUMPS, [LU](https://grokipedia.com/page/LU_decomposition), cuDSS) are selected and probed under policy and RAM constraints. 4. **Facade with lazy public API.** `import linstabpy` exposes analysis entry points without initializing [MPI](https://grokipedia.com/page/Message_Passing_Interface) or PETSc. 5. **Declarative configuration and BC recipes.** Case behavior comes from YAML and BC rules, not hardcoded geometry branches. 6. **Cost aware routing.** Anchor factorization, [preconditioner](https://grokipedia.com/page/Preconditioner) reuse, residual gates, break-even logic, and RAM autoselect turn single point resolvent into scalable campaigns. ### Technology stack overview | Layer | Choice | |-------|--------| | Language | Python ≥3.10, <3.14 | | Local arrays / sparse | NumPy, SciPy (CSC/CSR) | | Parallel runtime | MPI via `mpi4py` | | Distributed linear algebra | PETSc (`petsc4py`), **complex scalars required** | | Eigenproblems | SLEPc (`slepc4py`), EPS / weighted normal problems | | Sparse direct solvers | MUMPS (default via PCLU); optional cuDSS (GPU shell PC) | | Configuration | PyYAML; typed dataclasses + schema validation | | Packaging | conda-forge environment + editable pip install | | Outputs | NPZ result payloads; `metrics.json` + provenance | **Environment note:** `environment.yml` targets conda-forge **complex-scalar** PETSc/SLEPc. GPU PETSc is not a supported conda path; cuDSS is an optional, non-conda path. Plain `import linstabpy` does **not** initialize PETSc or MPI. ### High-level deployable mental model ```mermaid flowchart TB subgraph Entry A1["CLI: linstabpy.driver.cli"] A2["Library: import linstabpy"] end subgraph Pipeline["Staged HPC Pipeline"] direction TB P1["Config validate + typed policy"] P2["Baseflow load + preprocess"] P3["Mesh + SBP + filters"] P4["Coefficients + L0 + BC"] P5["Factorize shifted systems"] P6a["Stability EPS"] P6b["Resolvent H / gains"] P6c["Router beta-omega map"] P7["NPZ + metrics + provenance"] P1 --> P2 --> P3 --> P4 --> P5 P5 --> P6a P5 --> P6b P4 --> P6c P6c --> P6b P6a --> P7 P6b --> P7 P6c --> P7 end subgraph Backends B1["PETSc Mat/Vec/KSP"] B2["SLEPc EPS"] B3["MUMPS / LU / cuDSS"] end A1 --> P1 A2 --> P4 A2 --> P6a A2 --> P6b A2 --> P6c P5 --> B1 P5 --> B3 P6a --> B2 P6b --> B2 ``` ## System context Roles, scope, and external pieces are defined once in [overview](overview.md). The architecture-relevant cut is how the two entry points reach the runtime: ```mermaid flowchart LR subgraph Users U1["Researchers / Engineers"] U2["Library API consumers"] end subgraph Inputs YAML["YAML case config"] CFD["CFD baseflow npy/npz"] end subgraph LinStabPy CLI["driver.cli / run"] LIB["linstabpy public API"] end subgraph Runtime MPI["MPI / mpi4py"] PETSc["PETSc"] SLEPc["SLEPc"] MUMPS["MUMPS"] cuDSS["cuDSS optional"] NP["NumPy / SciPy"] end subgraph Outputs NPZ["Results NPZ"] MET["metrics.json + provenance"] end U1 --> CLI U2 --> LIB YAML --> CLI CFD --> CLI CLI --> MPI CLI --> PETSc CLI --> SLEPc PETSc --> MUMPS PETSc --> cuDSS CLI --> NP LIB --> PETSc CLI --> NPZ CLI --> MET ``` The CLI needs the full stack: MPI with PETSc and SLEPc plus MUMPS underneath, SciPy for local assembly, YAML and base-flow files in, NPZ results with metrics and provenance out. The library needs only NumPy and SciPy until a solve starts; PETSc enters lazily with the first distributed matrix or factorize call. Optional cuDSS factors ride as a PETSc preconditioner shell wherever the backend policy selects them. Solver vocabulary on this page is defined in the [solver glossary](solver-glossary.md). ## Container view Packages here are cohesive Python modules with a runtime and data contract. They run in one MPI process group as a single Python application. ### Domain module division | Domain | Type | Code paths | Role | |--------|------|------------|------| | **Case Orchestration** | Core | `driver/`, `__init__.py` | PETSc/MPI init, staged runs, CLI, payload, public lazy API | | **Configuration Management** | Supporting | `config/` | Merge/validate YAML·env·defaults; typed solver policy | | **Baseflow Management** | Core | `baseflow/` | Mean-flow ingest, sponge, shock cleanup, viscosity | | **Mesh Discretization** | Core | `mesh/` | Geometry, SBP ops, filters, deformation, masking | | **LNSE Operator Assembly** | Core | `operators/`, `boundary_conditions/` | Coefficients, L0/spatial assembly, BC recipes, Chu norm | | **Stability Analysis** | Core | `stability/`, `solvers/eigenvalue.py`, `spatial_eig.py` | Temporal/adjoint/spatial EPS + sensitivity | | **Resolvent Analysis** | Core | `resolvent/` | H operator, weighted gains, forcing models | | **Solver Routing** | Core | `solvers/router/` | \(\beta\)/\(\omega\) maps, anchors, residual gates, RAM autoselect | | **Linear Algebra** | Infrastructure | `linalg/`, `mumps_options.py` | Factor backends, shifted ops, cuDSS | | **Shared Utilities** | Infrastructure | `utils/`, `metrics_schema.py` | PETSc bridge, YAML, cost, provenance, Floquet | ### Domain module architecture ```mermaid flowchart TB subgraph Orchestration["Case Orchestration"] CLI["cli.py"] RUN["run.py"] PRE["preflight.py"] PAY["payload.py / io.py"] end subgraph Config["Configuration Management"] SCH["schema.py"] TYP["typed.py"] DEF["defaults.py"] end subgraph PhysicsPrep["Physics Preparation"] BF["baseflow"] MESH["mesh + SBP"] end subgraph Operators["LNSE Operator Assembly"] COEF["coefficients"] ASM["assembly"] SPAT["spatial"] BC["boundary_conditions"] EN["energy_norm"] end subgraph Analysis["Analysis Domains"] STAB["stability"] RES["resolvent"] RTR["solvers/router"] end subgraph Solvers["Solver Backends"] EIG["solvers/eigenvalue"] SEIG["solvers/spatial_eig"] end subgraph LA["Linear Algebra Infrastructure"] FAC["factorize"] BE["backends"] SH["shifted"] CU["cudss_factor"] end subgraph Utils["Shared Utilities"] PC["petsc_core / petsc_init"] YIO["yaml_io"] PROV["provenance / cost_metrics"] end CLI --> RUN RUN --> PRE RUN --> Config RUN --> BF RUN --> MESH RUN --> ASM RUN --> STAB RUN --> RES RUN --> RTR RUN --> PAY ASM --> COEF ASM --> BC ASM --> MESH ASM --> BF SPAT --> COEF STAB --> EIG STAB --> SEIG STAB --> FAC RES --> FAC RTR --> RES RTR --> FAC RTR --> ASM FAC --> BE FAC --> PC EIG --> FAC ``` ### Storage design LinStabPy is **file-oriented**, not database-backed. | Store | Form | Contract | |-------|------|----------| | Case policy | YAML on disk / env overrides | Merged into typed `Config`; schema-validated | | Mean flow | npy directory or npz | Primitives `RHO,U,V,W,T` on `(ny,nx)`; optional sponge | | Runtime operators | In-memory SciPy CSC or PETSc `Mat` | Sparse \(5 N_{gp} \times 5 N_{gp}\) (or reduced free DOFs) | | Results | NPZ | Fixed key contract in `driver.payload` | | Metrics | JSON | `metrics_schema.SCHEMA`: provenance, problem, cost, quality | | Provenance | Embedded in metrics / payload | Config hash/MD5, git SHA, hostname, factor backend/precision | There is no persistent application database. Reproducibility is achieved through config hashing, schema versioning, and immutable result contracts rather than transactional storage. ### Inter-domain module communication Communication is **in-process Python calls** (service calls, composition, and data hand-off of arrays/matrices), not network RPC. ```mermaid flowchart LR CO["Case Orchestration"] CM["Configuration"] BM["Baseflow"] MD["Mesh"] OA["LNSE Assembly + BC"] SA["Stability"] RA["Resolvent"] RR["Router"] LI["Linear Algebra"] UT["Utilities"] CO -->|"config dependency"| CM CO -->|"data"| BM CO -->|"service"| SA CO -->|"service"| RA CO -->|"service"| RR CO --> OA CO --> LI OA -->|"composition"| MD OA -->|"data"| BM OA --> LI SA -->|"service"| LI SA -->|"data"| OA RA -->|"service"| LI RA -->|"data"| OA RR -->|"service"| RA RR -->|"service"| LI RR -->|"data"| OA BM -->|"data"| MD MD --> UT SA --> UT RA --> UT LI --> UT ``` **Strength highlights** - Orchestration ↔ Stability/Resolvent: primary service calls - Assembly ↔ Mesh / Baseflow: composition and data (~strongest physics coupling) - Stability / Resolvent / Router ↔ Factorization: shared factorize-once spine **Boundary-conditions note:** The `boundary_conditions/` package is physically separate but **logically part of LNSE Operator Assembly**; BC application is the last non-interchangeable step of operator finalization. ## Component view ### Project structure ```text src/linstabpy/ ├── __init__.py # Lazy public API (no PETSc on import) ├── metrics_schema.py # Run metrics JSON contract ├── mumps_options.py # MUMPS option helpers ├── driver/ # CLI + case pipeline ├── config/ # Schema, typed Config, defaults ├── baseflow/ # BaseFlow + I/O + preprocess ├── mesh/ # Mesh, SBP, filters, FD weights ├── operators/ # LNSE assembly, coefficients, spatial, norms ├── boundary_conditions/ # Dirichlet/Neumann recipes ├── stability/ # Eigenvalue services, acceptance, sensitivity ├── resolvent/ # Operator, gain, analysis ├── solvers/ # EPS backends + router/ ├── linalg/ # factorize, backends, shifted, cuDSS └── utils/ # PETSc init/core, YAML, metrics, provenance ``` **Public library surface** (lazy exports from `linstabpy.__init__`): - `build_lnse_system`, `build_lnse_operator` - `compute_chu_norm` - `temporal_eigenvalues`, `adjoint_eigenvalues` - `wavemaker` - `resolvent_analysis`, `resolvent_map` - `BaseFlow` **CLI entry:** `linstabpy.driver.cli.main` → lazy import of `driver.run` (so `--version` / `--help` avoid PETSc). ### Core functional components #### Case orchestration (`driver`) | Component | Responsibility | |-----------|----------------| | `cli.py` | Argparse; cache check (config hash / git SHA); force-recompute env | | `run.py` | Full pipeline: init PETSc/MPI, load config/fields, mesh, assemble L0, stability/resolvent stages, payload | | `preflight.py` | Non-blocking semantic advisories (dissipation vs mesh coupling, freestream scale, etc.) | | `payload.py` / `io.py` | NPZ contract keys; save/load results | Notable runtime behaviors: - Stage wall-time and RSS delta logging; optional stage heartbeats - Env knobs: `LINSTABPY_CASE_CONFIG`, `LINSTABPY_CASE_OUTPUT_NPZ`, `LINSTABPY_FORCE_RECOMPUTE`, … - Alpha-filter treated as an **operator-level** property (not stage-dependent) #### Configuration (`config`) | Component | Responsibility | |-----------|----------------| | `schema.py` | Lenient allowed-key validation + FD-order checks → `ConfigValidationError` | | `typed.py` | Dataclasses: `Config`, `EigenCfg`, `LinSolveCfg`, `ResolventCfg`, `SpatialCfg`, `MeshCfg`, `MpiCfg` | | `defaults.py` | Dissipation and other apply defaults | **Resolution order:** YAML case → environment → defaults. #### Baseflow (`baseflow`) | Component | Responsibility | |-----------|----------------| | `BaseFlow` | `(ny, nx)` fields `RHO,U,V,W,T` + `Re,Ma,Pr,kappa,T_0`; viscosity via `sutherland_air` | | `grid_fields` | npy directory / npz load | | preprocess / sponge / shock_cleanup / viscosity | Readiness for linearization | #### Mesh (`mesh`) | Component | Responsibility | |-----------|----------------| | `Mesh` | Mutable runtime carrier; views `MeshGeometry`, `MeshOperators`, `MeshRuntimeState` | | SBP (`sbp.py`, `sbp_generator.py`, `fdweights.py`) | Differentiation, norms, dissipation, Mohseni–Colonius pole treatments | | `filters` | Mattsson-style stabilization attached to mesh runtime | PETSc-dependent mesh objects are lazy-imported. #### LNSE operator assembly (`operators` + `boundary_conditions`) **Canonical PETSc assembly order** (documented in `assembly.py`; **non-interchangeable**): 1. Assemble raw block operator from coefficients + FD stencils 2. Apply analytical `inv(B/1j)` left scaling → L0 3. Add sponge damping 4. Add artificial dissipation 5. **Last:** apply BC recipe (Dirichlet/Neumann); free-DOF sets State vector: **`q = [ρ, u, v, w, T]`**. Supports Cartesian and cylindrical coefficient paths; SciPy CSC and PETSc `Mat` backends. | Sub-module | Role | |------------|------| | Temporal Operator Assembly | Large sparse L0 for temporal/resolvent paths | | Physical Coefficients | LNSE coeffs, metrics, Chu energy-norm weights | | Spatial Operator Assembly | Spatial block operators reusing temporal helpers | | Boundary Condition Recipes | Declarative Dirichlet/Neumann rules and free-DOF masks | #### Stability (`stability` + `solvers`) | Sub-module | Role | |------------|------| | Stability Eigenvalue Services | `temporal_eigenvalues` / `adjoint_eigenvalues` / spatial entry; acceptance | | Eigen Solver Backends | SLEPc EPS shift-invert; residual reorder; conjugate-pair completeness | | Structural Sensitivity | `wavemaker` from direct + adjoint modes | #### Resolvent (`resolvent`) \[ H = C\,F\,(L_0 - (\gamma + i\omega) I)^{-1}\,F\,B \] - Factorize shifted operator once; reuse for forward/adjoint applies - Weighted normal operator shell for SLEPc; gains \(\sigma = \sqrt{\mathrm{Re}(\lambda)}\) - Optional conservative momentum forcing maps #### Solver router (`solvers/router`) | Piece | Role | |-------|------| | `resolvent_map` | Per-\(\beta\) L0 assembly; schedule frequency lines | | `band_marking` | Label STABLE / SMOOTH / NEAR_POLE bands | | `map_walk` / `anchors` | Real factor at anchors; later \(\omega\) uses factor as PC; residual gate → refactor | | `stable_sweep` | Stable-regime sweep strategies | | `recycling` / `break_even` | Cost-aware reuse decisions | | `ram_autoselect` | Backend choice under memory policy | | `certification` | Consistency / RSVD certification hooks | **Config nuance:** Router dissipation strength is read from top-level `dissipation_epsilon`, while non-router paths may use nested `operator_dissipation.epsilon`. Both enable flags and strength must be set or dissipation is a **silent no-op**. ### Technical support components #### Linear algebra (`linalg`) | Component | Responsibility | |-----------|----------------| | `factorize` / `Factorization` | solve, transpose solve, block solve, refactorize, destroy | | `backends` | Normalize/probe availability, precision, MUMPS BLR/ordering | | `shifted_operator` | Build \(L_0 - i\omega I\) and related shifts | | `cudss_factor` | Optional GPU factor as PETSc PC shell | **Default path:** KSPPREONLY + PCLU + MUMPS. #### Shared utilities (`utils` + `metrics_schema`) | Component | Responsibility | |-----------|----------------| | `petsc_init` / `petsc_core` | Init ordering (`mpi4py` before PETSc; thread level MULTIPLE when needed); SciPy↔PETSc conversion | | YAML I/O, cost metrics, provenance | Cross-cutting reproducibility and ops telemetry | | Floquet helpers | Supporting periodic-flow utilities | | `metrics_schema` | Versioned metrics JSON groups | ### Component responsibility division (summary) ```mermaid flowchart TB subgraph CoreBusiness ORCH["Orchestration: CLI run payload"] BF["Baseflow: ingest preprocess"] MESH["Mesh: SBP filters geometry"] OP["Operators+BC: L0 spatial norms"] ST["Stability: temporal adjoint spatial wavemaker"] RS["Resolvent: H gains modes"] RT["Router: beta-omega maps anchors"] end subgraph Infra CFG["Config: schema typed defaults"] LA["Linalg: factorize backends shifted"] UT["Utils: PETSc bridge metrics provenance"] end ORCH --> CFG ORCH --> BF ORCH --> MESH ORCH --> OP ORCH --> ST ORCH --> RS ORCH --> RT OP --> BF OP --> MESH ST --> OP ST --> LA RS --> OP RS --> LA RT --> OP RT --> RS RT --> LA LA --> UT MESH --> UT ST --> UT RS --> UT ``` ### Component interaction relationships - **Library consumers** may call assembly and analysis APIs directly, bypassing full CLI orchestration, but still depend on the same operator and factorization contracts. - **Router** is the only domain that systematically **re-assembles per \(\beta\)** and **walks many \(\omega\)** with measured cost feedback. - **Stability and Resolvent** never own mesh generation; they consume finalized operators and energy weights. - **Utilities** never own case flow; they underpin conversion, init, and telemetry. ## Key processes ### Core functional processes #### Temporal eigenvalue stability analysis (primary product path) End-to-end [hydrodynamic stability](https://grokipedia.com/page/Hydrodynamic_stability) for compressible flows: load mean flow, assemble LNSE + BCs, factorize shifted systems, compute temporal (and optionally adjoint) eigenvalues/eigenvectors. ```mermaid flowchart TD Start["Case Launch"] --> InitRuntime["Initialize PETSc and MPI"] InitRuntime --> ResolveConfig["Merge and Validate YAML Config"] ResolveConfig --> LoadBaseflow["Load CFD Mean-Flow Fields"] LoadBaseflow --> PreprocessBF["Preprocess Sponge Viscosity Shock Cleanup"] PreprocessBF --> BuildMesh["Build Structured Mesh and SBP Operators"] BuildMesh --> GenCoeffs["Compute LNSE Physical Coefficients"] GenCoeffs --> AssembleL0["Assemble Sparse L0 Operator"] AssembleL0 --> ApplyBCs["Apply Dirichlet Neumann BC Recipes"] ApplyBCs --> FactorShift["Factorize Shifted LNSE System"] FactorShift --> EigenSolve["SLEPc Temporal Adjoint Eigen Solve"] EigenSolve --> ValidateModes["Validate Modes Residuals Weights"] ValidateModes --> EmitPayload["Emit Spectra Modes Metrics Provenance"] EmitPayload --> End["Case Complete"] ``` ```mermaid sequenceDiagram participant CLI as driver_cli participant RUN as driver_run participant CFG as config participant BF as baseflow participant M as mesh participant OP as operators_BC participant LA as linalg participant ST as stability_solvers participant OUT as payload_io CLI->>RUN: run case RUN->>CFG: load validate typed Config RUN->>BF: load_case_fields preprocess viscosity RUN->>M: Mesh SBP filters RUN->>OP: coefficients L0 BC recipe RUN->>LA: factorize shifted L0 RUN->>ST: SLEPc temporal adjoint EPS ST->>LA: repeated solves RUN->>OUT: eigenvalues modes metrics provenance ``` **Key steps** 1. Runtime initialization, PETSc/MPI, logging, distributed orchestration 2. Configuration resolution, YAML/env/defaults; eigen/mesh/MPI policy 3. Baseflow ingest, npy/npz → `BaseFlow` + preprocess/transport 4. Mesh discretization, structured mesh, SBP (2, 4, 6, 8), filters, runtime state 5. Coefficient generation, Cartesian or cylindrical compressible LNSE coeffs 6. Operator assembly and BCs, [sparse](https://grokipedia.com/page/Sparse_matrix) L0 for \(q=[\rho,u,v,w,T]\); free DOFs 7. [Shift-invert](solver-glossary.md) factorization, MUMPS/LU/cuDSS factorize-once 8. Eigenanalysis, SLEPc [EPS](solver-glossary.md); residuals; optional structural sensitivity 9. Result emission, NPZ + metrics + provenance #### Resolvent gain analysis (single frequency) Shares prep through L0, then specializes in \((L_0 - i\omega I)\) factors and weighted gain extraction. ```mermaid flowchart TD Start["Resolvent Case Start"] --> InitRuntime["Initialize Runtime"] InitRuntime --> ResolveCfg["Resolve Frequency Weights Backend Policy"] ResolveCfg --> PrepMean["Provide Mesh-Aligned BaseFlow"] PrepMean --> SBPOps["Supply SBP Diff and Norm Operators"] SBPOps --> AssembleL0["Assemble L0 with BCs and Energy Maps"] AssembleL0 --> FactorOmega["Factor L0 minus i omega I"] FactorOmega --> BuildH["Build Resolvent H and H Hermitian Actions"] BuildH --> WeightedShell["Form Weighted Normal Operator Shell"] WeightedShell --> GainEPS["SLEPc Extract Gains and Mode Pairs"] GainEPS --> ResidualCheck["Validate Residuals"] ResidualCheck --> PackageOut["Package Gains Modes Metrics"] PackageOut --> End["Resolvent Complete"] ``` #### Beta–frequency resolvent gain map (cost-critical production path) ```mermaid flowchart TD Start["Map Campaign Launch"] --> ValidateMap["Validate Map Router Resource Policy"] ValidateMap --> LoopBeta["For Each Spanwise Beta"] LoopBeta --> AssembleBeta["Assemble LNSE Operator for Beta"] AssembleBeta --> MarkBands["Mark Frequency Bands"] MarkBands --> ScheduleWalk["Schedule Frequency Line Walk"] ScheduleWalk --> AnchorFactor["Factor Measured Anchors"] AnchorFactor --> WalkFreq["Walk Frequencies with PC Reuse"] WalkFreq --> ResidualGate{"Residual Gate OK?"} ResidualGate -->|Yes| ComputeGain["Compute Leading Resolvent Gain"] ResidualGate -->|No| Refactor["Refactor Anchor and Continue"] Refactor --> ComputeGain ComputeGain --> MoreFreq{"More Frequencies?"} MoreFreq -->|Yes| WalkFreq MoreFreq -->|No| MoreBeta{"More Betas?"} MoreBeta -->|Yes| LoopBeta MoreBeta -->|No| RecordMap["Record Map Cost Samples Certification"] RecordMap --> End["Map Complete"] ``` #### Spatial stability eigenproblem Config `solver.spatial` (fallback `solver.eigen`) → spatial block operators → generalized SLEPc EPS. ```mermaid flowchart TD Start["Spatial Eigen Start"] --> ResolveSpatial["Resolve solver.spatial Policy"] ResolveSpatial --> SupplyBF["Supply Mean-Flow Fields"] SupplyBF --> MeshIndex["Provide Mesh Indexing and SBP"] MeshIndex --> SpatialOps["Assemble Spatial Block Operators"] SpatialOps --> FactorBackend["Configure Shift-Invert Factor Backend"] FactorBackend --> SpatialEPS["Run Generalized SLEPc Spatial EPS"] SpatialEPS --> ReturnModes["Return Spatial Modes and Metadata"] ReturnModes --> End["Spatial Eigen Complete"] ``` #### Baseflow ingest and operator prep (shared foundation) ```mermaid flowchart TD Start["Prep Stage Start"] --> SchemaCheck["Validate Case Schema and FD Orders"] SchemaCheck --> LoadFields["Ingest npy or npz Primitives"] LoadFields --> Preprocess["Shock Cleanup Sponge Viscosity Fill"] Preprocess --> MeshCreate["Create or Deform Mesh Attach SBP Filters"] MeshCreate --> Coeffs["Generate LNSE Coefficients"] Coeffs --> BCRecipe["Apply BC Recipes Free-DOF Sets"] BCRecipe --> FinalizeL0["Finalize Sparse L0 for Downstream Solves"] FinalizeL0 --> End["Operator Ready"] ``` #### Adjoint and structural sensitivity (supporting path) ```mermaid flowchart TD Start["Direct Modes Available"] --> AdjointSolve["Run Adjoint Eigen Solve"] AdjointSolve --> CombineModes["Combine Direct and Adjoint Modes"] CombineModes --> Sensitivity["Compute Structural Sensitivity"] Sensitivity --> EmitSens["Emit Sensitivity Fields and Metadata"] EmitSens --> End["Sensitivity Complete"] ``` ### Technical processing workflows | Concern | Mechanism | |---------|-----------| | PETSc init ordering | `mpi4py` before PETSc; complex scalars mandatory for production | | Operator finalization | Fixed order: raw blocks → inv(B) scale → sponge → dissipation → BC last | | Shift construction | `linalg.shifted` builds \(L_0 - i\omega I\) (and related) once per target | | Eigen drive | SLEPc EPS shift-invert using factorization as PC/solve backend | | Resolvent drive | MatShell weighted normal operator; repeated H / Hᴴ applies via factors | | Map routing | Band mark → anchors → PC reuse → residual gate → optional stable sweep | | Backend choice | Probe + RAM autoselect among MUMPS/LU/cuDSS | ### Data flow paths ```text YAML + npy/npz → Config (typed) + BaseFlow (ny,nx fields) → Mesh (X,Y, Dx,Dy, W, filters) → Coefficients + Energy weights (Chu W_diag / invW_diag) → L0 PETSc Mat + free-DOF IS + BC-enforced rows → Factorization handle(s) → Spectra / gains / modes (distributed Vec → gathered arrays) → NPZ payload + metrics.json (provenance, cost, quality) ``` **Key data objects** | Object | Contract | |--------|----------| | `BaseFlow` | Parameters + `(ny,nx)` primitives; optional MU / dμ/dT | | `Mesh` / operators | `X,Y`, SBP `Dx,Dy,…`, weights `W`, optional dissipation | | `LNSESystem` / L0 | Sparse \(5 N_{gp} \times 5 N_{gp}\) (or reduced free DOFs) | | BC recipe | Rules → enforced rows + free-DOF index sets | | Energy weights | Chu norm diagonals `W_diag` / `invW_diag` | | Runtime NPZ | Fixed keys in `driver.payload` | | Metrics JSON | provenance, problem, cost, quality groups | ### Exception handling and quality gates | Gate | Behavior | |------|----------| | **Schema validation** | Unknown top-level keys / invalid FD orders → `ConfigValidationError` (fail fast before HPC cost) | | **Preflight advisories** | Non-blocking warnings (dissipation vs mesh coupling, freestream scale) | | **Residual gates (eigen)** | Mode residual checks, target distance, conjugate-pair completeness | | **Residual gates (resolvent/map)** | Force refactor when PC-reuse accuracy degrades | | **Acceptance services** | `stability/acceptance.py` post-checks before payload commit | | **Certification hooks** | Router consistency / RSVD-style certification metadata | | **Cache / recompute** | CLI config-hash and git-SHA cache checks; `LINSTABPY_FORCE_RECOMPUTE` override | There is no distributed saga or compensating transaction model: failures abort the stage with logged metrics where available; partial map results may still be recorded depending on campaign progress. ## Technical implementation ### Core module implementation #### Orchestration hub `driver/run.py` is a large (~1.8k line) façade that owns stage sequencing, timing, and payload emission. This concentration is intentional for scientific reproducibility (one place defines "a case") but is an architectural pressure point for change risk. Library users can bypass it by calling domain APIs directly. #### Operator assembly Assembly is the mathematical heart of the system: - Coefficients encode compressible LNSE physics (Cartesian/cylindrical). - SBP stencils supply discrete derivatives and norms with summation-by-parts stability properties. - Left scaling by analytical `inv(B/1j)` produces the [temporal operator](core-concepts.md) L0 used by eigen and resolvent paths. - Sponge and artificial dissipation are operator-level stabilizations. - BC recipes mutate rows/columns and publish free-DOF sets **last**, so solvers never see unconstrained boundary rows. #### Factorization layer `Factorization` is the performance spine: - Default: **KSPPREONLY + PCLU + MUMPS** - API: solve, transpose solve, block solve, refactorize, destroy - Optional cuDSS as a PETSc PC shell for GPU-capable deployments - Shifted systems built once per \((\beta,\omega)\) (or eigen shift) and reused across SLEPc iterations and resolvent matvecs #### Analysis products - **Stability:** SLEPc EPS on shifted L0; adjoint path; wavemaker from direct×adjoint - **Resolvent:** \(H = C F (L_0 - (\gamma + i\omega) I)^{-1} F B\); weighted normal eigenproblem; \(\sigma = \sqrt{\mathrm{Re}(\lambda)}\) - **Router:** per-\(\beta\) assembly + frequency line walks with measured anchors and residual-gated reuse ### Key algorithm design | Algorithm area | Design choice | Rationale | |----------------|---------------|-----------| | Spatial discretization | SBP FD orders 2, 4, 6, 8 + Mattsson-style filters | Energy-stable structured discretization; MATLAB LinStab2D parity | | Axisymmetric poles | Mohseni–Colonius treatments | Correct pole behavior on cylindrical grids | | Temporal eigen | Shift-invert via direct factors | Dense spectral regions near targets; robust on stiff LNSE | | Resolvent gains | Weighted normal MatShell + SLEPc | Avoid forming dense H; reuse factors for H and Hᴴ | | Map campaigns | Anchor factors as PC + residual gate | Amortize factorization across nearby \(\omega\); protect accuracy | | Backend policy | RAM autoselect + MUMPS BLR options | Fit large free-DOF systems into node memory | | Energy inner product | Chu norm diagonals | Physically meaningful compressible disturbance energy | ### Data structure design | Structure | Representation | Notes | |-----------|----------------|-------| | Mean-flow fields | NumPy `(ny, nx)` arrays on `BaseFlow` | Mesh-aligned primitives | | Mesh operators | Sparse differentiation/norm matrices | Attached to `Mesh` runtime views | | L0 / spatial blocks | SciPy CSC **or** PETSc `Mat` | Production path is PETSc distributed | | Free DOFs | PETSc index sets / masks | After BC enforcement | | Factors | Opaque `Factorization` handle | Owns KSP/PC lifecycle | | Modes / gains | Distributed `Vec` → NumPy in payload | Gathered for NPZ contract | | Map records | Structured arrays + cost samples | Router recording + certification | State dimension scales as **5 degrees of freedom per grid point** (density, three velocities, temperature), reduced by Dirichlet eliminations. ### Performance optimization strategies 1. **Factorize-once, solve-many**, Dominant cost is sparse direct factorization; eigen iterations and resolvent applies reuse factors. 2. **Anchor + PC reuse on maps**, Pay full factor cost at anchors; nearby frequencies reuse as right preconditioners until residual gate trips. 3. **Break-even / recycling logic**, Router estimates when refactor is cheaper than struggling iterative reuse. 4. **RAM autoselect**, Choose MUMPS/LU/cuDSS under memory policy; MUMPS BLR compression where available. 5. **Lazy imports**, Keep CLI help and library import off the PETSc critical path. 6. **Stage metrics**, Wall time and RSS deltas identify which stage dominates (assembly vs factor vs EPS). 7. **Optional parallelism across \(\beta\)**, Conceptual extension: parallelize independent \(\beta\) lines when memory allows; keep per-line \(\omega\) walks sequential for factor reuse. 8. **SciPy local path**, Useful for small debugging assemblies without full MPI cost. **Primary bottlenecks** - Sparse direct factor memory and time (grows superlinearly with free DOFs) - Per-\(\beta\) reassembly on large maps - Monolithic `run.py` and large assembly/coefficient modules (developer throughput, not FLOPs) - Dual dissipation config paths (silent no-op risk → wasted or wrong HPC runs) ### Scalability design | Axis | Strategy | |------|----------| | **Grid size** | MPI-distributed PETSc matrices; direct sparse solvers | | **Parameter space (\(\beta,\omega\))** | Router map campaigns with cost-aware reuse | | **Analysis modes** | Same L0 family → many products without separate physics stacks | | **Hardware** | CPU complex PETSc/MUMPS default; optional cuDSS GPU factors | | **Embedding** | Library API for pipeline integration; CLI for batch cases | | **Extension points** | New BC recipes; new factor backends via registry; new router band strategies; spatial/temporal coefficient models | ### Security design LinStabPy is a **trusted scientific HPC toolkit**, not a multi-tenant service. Security posture is appropriate to that model: | Concern | Approach | |---------|----------| | Multi-tenant auth / isolation | **Out of scope**, runs under cluster user accounts and MPI job schedulers | | Input trust | YAML and npy/npz are trusted researcher inputs; schema validation rejects unknown keys and bad FD orders but is not a sandbox | | Code execution | No dynamic code execution from case files beyond configured numeric policy | | Secrets | No cloud credentials subsystem; provenance records host/git/config hash, not secrets | | Supply chain | conda-forge PETSc/SLEPc pins; editable package install in controlled envs | | Result integrity | Config hash, git SHA, metrics schema, residual/certification gates support scientific integrity more than adversarial security | Hardening guidance for operators: run inside standard HPC job isolation; treat case directories as trusted; do not expose the CLI as an unauthenticated network service. ## Deployment architecture ### Runtime environment requirements | Requirement | Detail | |-------------|--------| | Python | ≥3.10, <3.14 | | MPI | Working MPI + `mpi4py` | | PETSc/SLEPc | **Complex scalar** builds via petsc4py/slepc4py | | MUMPS | Default sparse-direct backend through PETSc | | NumPy/SciPy/PyYAML | Local numerics and config | | Optional cuDSS | GPU sparse direct; not the supported conda default path | | Resources | Memory-bound for large free-DOF factors; multi-rank for distributed Mat | **Init constraints** - Initialize `mpi4py` before PETSc - Prefer thread level MULTIPLE when backends require it - Do not expect `import linstabpy` alone to set up the parallel runtime ### Deployment topology structure ```mermaid flowchart TB subgraph UserSpace["User / CI workspace"] YAML["Case YAML"] DATA["Baseflow npy or npz"] OUT["NPZ + metrics.json"] end subgraph Launch["Launch layer"] SCHED["HPC scheduler mpirun sbatch"] CLI["linstabpy.driver.cli"] end subgraph Ranks["MPI process group"] R0["Rank 0 orchestration logging"] R1["Rank 1 PETSc ownership"] Rn["Rank N distributed Mat Vec"] end subgraph Native["Native solver stack"] PETSc["PETSc"] SLEPc["SLEPc"] MUMPS["MUMPS"] GPU["Optional cuDSS"] end YAML --> CLI DATA --> CLI SCHED --> CLI CLI --> R0 CLI --> R1 CLI --> Rn R0 --> PETSc R1 --> PETSc Rn --> PETSc PETSc --> MUMPS PETSc --> GPU PETSc --> SLEPc R0 --> OUT ``` **Typical deployment unit:** one MPI job = one case (or one map campaign). There is no long-running application server. **Library embedding topology:** parent Python process imports `linstabpy`, explicitly initializes PETSc/MPI when needed, and calls public analysis APIs; I/O may be owned by the parent pipeline. ### Scalability in deployment - **Strong scaling:** PETSc distribution of sparse operators and vectors across ranks - **Parameter scaling:** router \(\beta\) loops and frequency walks; future multi-\(\beta\) parallelism where memory permits - **Vertical scaling:** larger-memory nodes for MUMPS factors; BLR compression; optional GPU factors - **Throughput scaling:** many independent cases via scheduler array jobs (embarrassingly parallel across cases) ### Monitoring and operations | Signal | Source | Use | |--------|--------|-----| | Stage wall time | Driver stage logging | Find slow prep vs factor vs EPS | | RSS / peak memory | Cost metrics | Backend and mesh sizing | | Optional GPU memory | Metrics when enabled | cuDSS capacity planning | | Residuals / acceptance | Stability & resolvent validators | Scientific quality gate | | Map cost samples | Router recording | Tune anchors and break-even | | Provenance | Config hash, git SHA, hostname, backend/precision | Reproduce and audit runs | | Preflight advisories | `preflight.py` | Catch policy/mesh mismatches early | | Heartbeats | Optional stage heartbeats | Long-running map health on HPC | **Operations guidance** 1. Validate YAML with schema before large allocations. 2. Start with single-frequency or single-shift cases before full \(\beta\)–\(\omega\) maps. 3. Watch factor memory first; adjust ranks, BLR, or backend via RAM policy. 4. Ensure dissipation enable + strength are set consistently (router top-level vs nested operator paths). 5. Archive NPZ + `metrics.json` together for provenance. 6. Use config-hash caching and `LINSTABPY_FORCE_RECOMPUTE` deliberately in CI or parameter sweeps. ### Development guidance | Audience | Guidance | |----------|----------| | **New contributors** | Learn the stage order and assembly order first; treat BC-last and factorize-once as invariants | | **Physics developers** | Extend `operators/coefficients` and BC recipes; keep SciPy and PETSc paths aligned | | **Solver developers** | Extend `linalg/backends` and router walk strategies; never bypass residual gates silently | | **API consumers** | Prefer lazy public exports; initialize PETSc only on real solve paths | | **Maintainers** | Consider splitting `run.py` into stage modules behind one façade; unify dissipation config; publish short ADRs for assembly order, factor policy, and router gates | ### Architecture risks and recommendations | Item | Severity | Notes / recommendation | |------|----------|------------------------| | No checked-in ADRs | Medium | Publish ADRs for assembly ordering, factor backend policy, residual gates, dissipation single source of truth | | Monolithic `driver/run.py` | Medium | Extract prep / assemble / stability / resolvent / emit stages; keep one façade | | Large assembly / coefficients / router modules | Medium | Add internal seams without breaking public contracts | | Dual dissipation config paths | Medium | Unify nested fields; keep aliases with deprecation warnings | | Typed config still migrating | Info | Consolidate `typed.py` with `utils._solver_cfg` dual resolution | | Spatial path less CLI-central | Low | Document spatial and map entry points beside lazy public API | | GPU path optional / non-conda | Info | Keep CPU complex PETSc as supported default; document cuDSS clearly | ## Summary LinStabPy implements a **clean, domain-aligned LNSE analysis architecture**: - YAML-driven **Case Orchestration** feeds mesh-aligned **Baseflow** into **SBP Mesh** discretization. - **LNSE Operator Assembly** (with declarative **BC recipes**) produces a single family of BC-enforced operators. - A shared **Linear Algebra** factorization layer underpins **Stability**, **Resolvent**, and cost-aware **Router** campaigns on **PETSc/SLEPc**. - Package boundaries match scientific responsibilities; external CFD generation, GUIs, and unstructured 3D frameworks remain firmly out of scope. The main architectural pressures are **orchestration monolith size**, **dual config paths for dissipation/solver policy**, and **missing formal ADRs**, not a mismatch between intended domains and the codebase layout. --- *LinStabPy, linear stability and resolvent analysis of compressible flows.* --- ## Workflow Source: https://linstabpy.vercel.app/workflow/ # Core workflows !!! info "Conventions" $\omega = \mathrm{eig}(i L_0)$, $\exp(-i\omega t)$ in time, unstable exactly when $\mathrm{Im}(\omega) > 0$. [Resolvent](https://grokipedia.com/page/Resolvent_formalism) $H = C F (L_0 - (\gamma + i\omega) I)^{-1} F B$. [SBP](https://grokipedia.com/page/Summation_by_parts) orders $2$, $4$, $6$, $8$. Default shift solve `KSPPREONLY` with `PCLU` through [MUMPS](https://en.wikipedia.org/wiki/MUMPS_(software)). Assembly order is fixed with BC rows last. ## Workflow overview ### System positioning LinStabPy is a Python CLI/library toolkit for **linear stability and resolvent analysis of [compressible](https://grokipedia.com/page/Compressible_flow) flows**. It ports LinStab2D-style LNSE workflows to [MPI](https://grokipedia.com/page/Message_Passing_Interface)-enabled HPC. Operationally it is a **staged scientific pipeline**: **config → baseflow → mesh/SBP → coefficients → L0 + BCs → factorize → eigen/resolvent → payload** The goal is to turn mean flow dumps into stability modes, resolvent [gains](core-concepts.md), and sensitivity results for [transition](http://www.scholarpedia.org/article/Transition_to_turbulence) prediction, receptivity, and flow-control design, without MATLAB, on reproducible YAML-driven MPI runs. **In scope:** baseflow ingest/preprocess, structured SBP meshes, LNSE assembly and BC recipes, temporal/spatial/[adjoint](https://grokipedia.com/page/Adjoint) eigenproblems, resolvent gains and β–ω maps, [sparse](https://grokipedia.com/page/Sparse_matrix) factorization (MUMPS/[LU](https://grokipedia.com/page/LU_decomposition)/optional [cuDSS](https://docs.nvidia.com/cuda/cudss/)), cost-aware solver routing, CLI orchestration, provenance. **Out of scope:** nonlinear [CFD](https://grokipedia.com/page/Computational_fluid_dynamics) that produces the mean flow, GUIs/visualization products, unstructured or full 3D volume frameworks, multiphysics outside LNSE, cloud multi-tenancy. ### System main workflows | Workflow | Role | Importance | Primary entry | |---|---|---|---| | **Temporal eigenvalue stability analysis** | Primary scientific product: spectra, modes, optional adjoint/sensitivity | Highest | `driver/run.py`, `stability._eigenvalues` | | **Resolvent gain analysis** | Optimal forcing/response and energy gains at target ω | Very high | `resolvent.analysis`, driver stages | | **Beta–frequency resolvent gain map** | Production parametric HPC campaign over (β, ω) | Very high (cost-critical) | `solvers/router/resolvent_map.py` | | **Spatial Stability Eigenproblem** | Streamwise wavenumber / spatial growth | High | `solvers/spatial_eig.py`, spatial assembly | | **Baseflow ingest and operator prep** | Shared foundation for all analysis products | High (shared spine) | driver preflight/setup | | **Adjoint and structural sensitivity** | Control/receptivity extension of temporal path | Supporting | `stability` adjoint + `wavemaker` | ### Core execution paths There are two complementary execution surfaces: 1. **CLI / case driver path** (`src/linstabpy/driver/run.py`), heavy runtime hub: [PETSc](https://en.wikipedia.org/wiki/PETSc)/MPI init, staged logging/heartbeats, config resolution, distributed vectors, end-to-end orchestration, payload emission. 2. **Library API path** (`import linstabpy`), lightweight: **does not** start PETSc/MPI on import; public symbols load lazily (`temporal_eigenvalues`, `resolvent_analysis`, `resolvent_map`, `build_lnse_system`, `BaseFlow`, etc.). Both paths converge on the same operator core and factorization spine. ```mermaid flowchart TD EntryCLI["CLI or driver run.py"] EntryLib["Library lazy API"] SharedPrep["Shared prep: config baseflow mesh L0 BCs"] BranchEigen["Temporal spatial adjoint eigen"] BranchRes["Single-point resolvent"] BranchMap["Beta frequency gain map"] FactorCore["Factorize-once solve-many linalg"] PayloadOut["Payload metrics provenance"] EntryCLI --> SharedPrep EntryLib --> SharedPrep SharedPrep --> BranchEigen SharedPrep --> BranchRes SharedPrep --> BranchMap BranchEigen --> FactorCore BranchRes --> FactorCore BranchMap --> FactorCore FactorCore --> PayloadOut ``` ### Key process nodes | Node | Domain | Contract | |---|---|---| | Runtime init | Case Orchestration | PETSc/MPI, rank-0 logging, stage timers, heartbeats | | Config merge/validate | Configuration | YAML + env + defaults → typed policy; schema gate | | BaseFlow load/preprocess | Baseflow | npy/npz → mesh-aligned primitives + transport | | Mesh + SBP | Mesh | Geometry, even-order SBP (2, 4, 6, 8), filters, masking, poles | | Coefficients | LNSE Assembly | Cartesian/cylindrical LNSE coeffs from mean flow | | L0 + BC recipes | LNSE Assembly | Sparse operator for `q = [ρ, u, v, w, T]`; free-DOF sets | | Factorization | Linear Algebra | Default KSPPREONLY + PCLU + MUMPS; optional cuDSS/iterative | | SLEPc EPS / gain | Stability / Resolvent | Shift-invert eigen or weighted normal operator | | Router walk | Solver Routing | Anchors, PC reuse, residual gates, certification | | Payload | Case Orchestration | Spectra/modes/gains, cost metrics, provenance NPZ | ### Process coordination mechanisms - **Single orchestration hub:** `driver/run.py` sequences stages; domains do not call each other ad hoc across the full case. - **One operator family, many products:** temporal, adjoint, spatial, resolvent, and maps all consume BC-enforced LNSE operators (spatial reuses temporal prep helpers). - **Factorize-once, solve-many:** shared performance spine for [shift-invert](solver-glossary.md) eigen and resolvent \(H\) actions. - **Config as reproducibility gate:** invalid structure/FD orders fail before expensive factors. - **MPI coordination:** distributed Mat/Vec via PETSc; rank-0 owns human-readable logs and many I/O decisions. - **Router feedback loop:** measured factor time/memory/compression samples inform reuse, backend choice, and certification on map campaigns. ## Main workflows ### Baseflow ingest and operator prep (shared foundation) **Why it matters:** Converts on-disk CFD mean flow into a mathematically consistent, BC-ready discrete LNSE operator. Every analysis product depends on this path. **Entry:** driver preflight/setup stages; library callers that build systems via `build_lnse_system` / related assembly APIs. #### Execution order and dependencies ```mermaid flowchart TD PrepStart["Prep stage start"] SchemaCheck["Validate case schema and FD orders"] LoadFields["Ingest npy or npz primitives"] Preprocess["Shock cleanup sponge viscosity fill"] MeshCreate["Create or deform mesh attach SBP filters"] Coeffs["Generate LNSE coefficients"] BCRecipe["Apply BC recipes free-DOF sets"] FinalizeL0["Finalize sparse L0 SciPy or PETSc"] PrepEnd["Operator ready"] PrepStart --> SchemaCheck SchemaCheck --> LoadFields LoadFields --> Preprocess Preprocess --> MeshCreate MeshCreate --> Coeffs Coeffs --> BCRecipe BCRecipe --> FinalizeL0 FinalizeL0 --> PrepEnd ``` | Step | Operation | Inputs | Outputs | Code anchors | |---|---|---|---|---| | 1 | Schema validation | Merged case dict | Pass or `ConfigValidationError` | `config/schema.py` | | 2 | Field I/O | npy dir or npz | Dict of 2-D float fields | `baseflow/grid_fields.py` | | 3 | Preprocess/transport | Raw primitives | Cleaned BaseFlow + μ, sponge | `preprocess`, `sponge`, `viscosity` | | 4 | Mesh/SBP | x, y, FD order, deformation | `Mesh` + operators + runtime state | `mesh/mesh.py`, `mesh/sbp.py` | | 5 | Coefficients | BaseFlow + metrics | Physical LNSE coefficient blocks | `operators/coefficients.py` | | 6 | BC recipes | YAML-style rules, mesh indices | Enforced rows/cols, free-DOF sets | `boundary_conditions/recipe.py`, `dirichlet.py` | | 7 | Finalize L0 | Coeffs, FD mats, sponge, backend request | SciPy CSC or PETSc `Mat` L0 | `operators/assembly.py` | **State transitions:** `raw files → validated config → BaseFlow (ny,nx) → MeshRuntimeState → coefficient tensors → sparse L0 (constrained) → free-DOF mask`. **Notes for implementers:** - BaseFlow fields are contractually **(ny, nx)** on the computational mesh. - SBP covers even **orders 2, 4, 6, 8**, including Mohseni–Colonius pole treatments for axisymmetric grids. - Boundary conditions are a separate package but **functionally part of operator assembly**. - SciPy CSC is typical for local assembly; PETSc Mat is used for distributed solve paths (`utils/petsc_core.py` bridge). ### Temporal eigenvalue stability analysis (primary product path) **Why it matters:** [Hydrodynamic stability](https://grokipedia.com/page/Hydrodynamic_stability) spectra and [eigenvectors](https://grokipedia.com/page/Eigenvalues_and_eigenvectors) for transition and control studies; reference orchestration pattern for other modes. **Entry:** `src/linstabpy/driver/run.py` or `temporal_eigenvalues` / stability APIs. #### End-to-end flow ```mermaid flowchart TD StartNode["Case launch"] InitRuntime["Initialize PETSc and MPI"] ResolveConfig["Merge and validate YAML config"] LoadBaseflow["Load CFD mean-flow fields"] PreprocessBF["Preprocess sponge viscosity shock cleanup"] BuildMesh["Build structured mesh and SBP operators"] GenCoeffs["Compute LNSE physical coefficients"] AssembleL0["Assemble sparse L0 operator"] ApplyBCs["Apply Dirichlet Neumann BC recipes"] FactorShift["Factorize shifted LNSE system"] EigenSolve["SLEPc temporal adjoint eigen solve"] ValidateModes["Validate modes residuals weights"] EmitPayload["Emit spectra modes metrics provenance"] EndNode["Case complete"] StartNode --> InitRuntime InitRuntime --> ResolveConfig ResolveConfig --> LoadBaseflow LoadBaseflow --> PreprocessBF PreprocessBF --> BuildMesh BuildMesh --> GenCoeffs GenCoeffs --> AssembleL0 AssembleL0 --> ApplyBCs ApplyBCs --> FactorShift FactorShift --> EigenSolve EigenSolve --> ValidateModes ValidateModes --> EmitPayload EmitPayload --> EndNode ``` #### Key technical process details 1. **Runtime initialization** - `init_petsc_from_argv()` then MPI/PETSc imports in the driver. - Stage wall-clock and RSS deltas (`_stage_start` / `_log_stage`); optional heartbeats via `LINSTABPY_STAGE_HEARTBEAT_SECONDS` (default 30s). - Config path from `LINSTABPY_CASE_CONFIG` (and legacy aliases) or `configs/parity_config.yaml`. 2. **Configuration resolution** - Merge YAML, environment, defaults into typed objects (`config/typed.py`). - Schema rejects unknown top-level structure and invalid FD orders while remaining lenient for legacy flat keys. 3. **Shared prep** - Same baseflow → mesh → coefficients → L0 + BC path as in [Baseflow ingest and operator prep](#baseflow-ingest-and-operator-prep-shared-foundation). - Driver uses `load_case_fields`, `Mesh`, `build_lnse_system_result`, BC analysis helpers. 4. **Shift-invert factorization** - Factorization layer (`linalg/factorize.py`): **direct-by-default** factorize-once, solve-many. - Default: **KSPPREONLY + PCLU + MUMPS**. - Iterative policy is **explicit** (not silent fallback) for memory-bound cases. - Optional cuDSS shelled as PETSc PC when GPU factors are selected. - Interfaces: solve, transpose solve, block solve, refactorize, destroy. 5. **Eigenanalysis** - [SLEPc](https://en.wikipedia.org/wiki/SLEPc) [EPS](solver-glossary.md) via `stability/_eigenvalues.py` and `solvers/eigenvalue.py`. - Temporal and adjoint paths; residual reordering; safe PETSc info access. - Acceptance helpers: target-distance postcheck, conjugate-pair completeness, stability acceptance (`stability/acceptance.py`). 6. **Result emission** - `build_runtime_payload` packages eigenvalues, eigenvectors, metrics, provenance. - Output NPZ path from env or `{prefix}_{config_md5}_{git_sha8}.npz`. #### Input / output data flow | Direction | Artifacts | |---|---| | **In** | YAML case; npy/npz mean flow (`ρ,u,v,w,T`, mesh, sponge); MPI ranks; factor backend policy | | **Internal** | BaseFlow, Mesh/SBP, L0, free-DOF sets, factored shifted operators, EPS workspace | | **Out** | Eigenvalues/eigenvectors, residual metrics, acceptance flags, cost/RSS stage logs, provenance (config hash, git SHA) | #### Dependencies ``` Configuration → Baseflow → Mesh → Operator Assembly → Linear Algebra → Stability → Payload ``` Orchestration **service-calls** Stability and **module-depends** on Operator Assembly and Linear Algebra. ### Resolvent gain analysis **Why it matters:** Amplifier/receptivity analysis via optimal forcing and response at chosen frequencies. **Mathematical core:** \[ H = C\,F\,(L_0 - (\gamma + i\omega) I)^{-1}\,F\,B \] Gains from the weighted eigenproblem on forms such as \(H^{H} W H\) (energy weights \(W\), often [Chu](core-concepts.md) norm), implemented as PETSc **[MatShell](solver-glossary.md)** + SLEPc. #### Flow ```mermaid flowchart TD ResStart["Resolvent case start"] InitRT["Initialize runtime"] ResCfg["Resolve frequency weights backend policy"] PrepMean["Provide mesh-aligned BaseFlow"] SBPOps["Supply SBP diff and norm operators"] AsmL0["Assemble L0 with BCs and energy maps"] FactorOm["Factor L0 minus i omega I"] BuildH["Build resolvent H and H Hermitian actions"] WShell["Form weighted normal operator shell"] GainEPS["SLEPc extract gains and mode pairs"] ResCheck["Validate residuals"] PkgOut["Package gains modes metrics"] ResEnd["Resolvent complete"] ResStart --> InitRT InitRT --> ResCfg ResCfg --> PrepMean PrepMean --> SBPOps SBPOps --> AsmL0 AsmL0 --> FactorOm FactorOm --> BuildH BuildH --> WShell WShell --> GainEPS GainEPS --> ResCheck ResCheck --> PkgOut PkgOut --> ResEnd ``` #### Key nodes | Step | Detail | Code | |---|---|---| | Policy | ω, weights, backend, NCV-related sizing | `config/typed.py`, resolvent sections | | L0 + energy | Same assembly spine; `compute_chu_norm` for compressible energy inner product | `operators/assembly.py`, `energy_norm.py` | | Factor (L0 − iωI) | Reused across many H applies | `linalg/factorize.py` | | H / Hᴴ | Input map B, output map C; PETSc-aware vectors | `resolvent/operator.py` | | Weighted normal | `invW * Hᴴ * W * H(v)` MatShell | `resolvent/gain.py` | | Analysis driver | Optimal pairs, gains, residual norms; MATLAB-aligned path | `resolvent/analysis.py` | **I/O:** Inputs extend temporal prep with frequency list/weights and forcing/response maps; outputs are gains, forcing/[response mode](core-concepts.md) shapes, residual certificates, solver metrics. ### Beta–frequency resolvent gain map (cost-critical campaign) **Why it matters:** Production-scale parametric studies of leading resolvent gain over spanwise [wavenumber](https://grokipedia.com/page/Wavenumber) β and frequency ω, with **cost-aware** factorization reuse. **Entry:** `src/linstabpy/solvers/router/resolvent_map.py` (public `resolvent_map`). #### Campaign flow ```mermaid flowchart TD MapStart["Map campaign launch"] ValMap["Validate map router resource policy"] LoopBeta["For each spanwise beta"] AsmBeta["Assemble LNSE operator for beta"] MarkBands["Mark frequency bands"] SchedWalk["Schedule frequency line walk"] AnchorFac["Factor measured anchors"] WalkFreq["Walk frequencies with PC reuse"] ResGate{"Residual gate OK?"} CompGain["Compute leading resolvent gain"] Refactor["Refactor anchor and continue"] MoreFreq{"More frequencies?"} MoreBeta{"More betas?"} RecMap["Record map cost samples certification"] MapEnd["Map complete"] MapStart --> ValMap ValMap --> LoopBeta LoopBeta --> AsmBeta AsmBeta --> MarkBands MarkBands --> SchedWalk SchedWalk --> AnchorFac AnchorFac --> WalkFreq WalkFreq --> ResGate ResGate -->|Yes| CompGain ResGate -->|No| Refactor Refactor --> CompGain CompGain --> MoreFreq MoreFreq -->|Yes| WalkFreq MoreFreq -->|No| MoreBeta MoreBeta -->|Yes| LoopBeta MoreBeta -->|No| RecMap RecMap --> MapEnd ``` #### Routing mechanics (implementation-aligned) 1. **Per-β operator assembly**, one BC-enforced LNSE operator per spanwise wavenumber (β enters the linearized physics). 2. **Band marking**, labels such as NEAR_POLE, SMOOTH, STABLE (`band_marking.py`) guide walk strategy. 3. **Line walk (`map_walk.LineWalker`)** - First frequencies on a line pay a **real factor** and measure cost. - Later points apply \(K(\omega)^{-1}\) using that factor as a **right [preconditioner](https://grokipedia.com/page/Preconditioner)**. - **Residual gate** forces a new factor when accuracy degrades; walks may abort or defer (`WalkDeferred`). 4. **Stable-regime sweeps (`stable_sweep.py`)**, anchor-preconditioned sweeps across nearby frequencies; records factor time, memory, compression for routing feedback. 5. **Backend policy**, `linalg/backends.py` normalizes/probes MUMPS/LU/cuDSS; RAM autoselect and break-even logic influence choices. 6. **Leading gain**, still from SLEPc/resolvent analysis at each scheduled point. 7. **Recording**, map arrays, cost keys (wall assemble/factor/solve, peak RSS/GPU, ranks, git SHA, etc.), certification metadata; atomic save patterns in `map_record`. #### Dissipation coordination (map-specific operational rule) Whether SBP sawtooth dissipation is needed depends on the **mesh**: - On **rectilinear** tensor meshes, metric cross-terms vanish; dissipation can be a costly no-op on gain accuracy. - On **curvilinear-in-both-directions** meshes, undamped operators can show spurious near-singular resonances. - Router expects **both**: `numerics.operator_dissipation.enabled` (build Q_diss on mesh) **and** top-level `dissipation_epsilon` (strength subtracted in `build_lnse_system`). Enabling only one is a **silent no-op**. #### Sequence: orchestration vs router vs linear algebra ```mermaid sequenceDiagram participant Driver as CaseOrchestration participant Cfg as Configuration participant Asm as LNSEAssembly participant Router as ResolventMapRouter participant Walk as LineWalker participant LA as FactorizationLayer participant Res as ResolventAnalysis participant Out as PayloadRecord Driver->>Cfg: Validate map router resource policy Driver->>Router: Launch resolvent_map campaign loop Each beta Router->>Asm: Assemble L0 for beta plus BCs Router->>Router: Mark frequency bands Router->>Walk: Schedule line walk Walk->>LA: Factor anchor measure cost loop Frequencies on line Walk->>LA: Solve with factor as PC or refactor Walk->>Res: Leading gain via SLEPc Walk->>Walk: Residual gate decide reuse end end Router->>Out: Map gains cost samples certification ``` ### Spatial stability eigenproblem **Why it matters:** Spatial growth/decay and streamwise wavenumber analysis via generalized eigenproblems. ```mermaid flowchart TD SpStart["Spatial eigen start"] SpCfg["Resolve solver.spatial policy"] SpBF["Supply mean-flow fields"] SpMesh["Provide mesh indexing and SBP"] SpOps["Assemble spatial block operators"] SpFac["Configure shift-invert factor backend"] SpEPS["Run generalized SLEPc spatial EPS"] SpOut["Return spatial modes and metadata"] SpEnd["Spatial eigen complete"] SpStart --> SpCfg SpCfg --> SpBF SpBF --> SpMesh SpMesh --> SpOps SpOps --> SpFac SpFac --> SpEPS SpEPS --> SpOut SpOut --> SpEnd ``` | Step | Behavior | |---|---| | Config | `solver.spatial` with **fallback to `solver.eigen`** when spatial YAML absent | | Assembly | `operators/spatial.py` reuses temporal preparation/index helpers; SciPy and PETSc block-matrix paths | | Solve | `solvers/spatial_eig.py`, generalized SLEPc EPS, direct shift-invert backend selection akin to temporal | | Output | Spatial eigenvalues/modes and solver metadata | ### Adjoint and structural sensitivity (supporting path) Extends temporal eigenanalysis without a separate full case pipeline: ```mermaid flowchart TD SensStart["Direct modes available"] AdjSolve["Run adjoint eigen solve"] Combine["Combine direct and adjoint modes"] SensComp["Compute structural sensitivity wavemaker"] SensEmit["Emit sensitivity fields and metadata"] SensEnd["Sensitivity complete"] SensStart --> AdjSolve AdjSolve --> Combine Combine --> SensComp SensComp --> SensEmit SensEmit --> SensEnd ``` - Adjoint path: `adjoint_eigenvalues` (lazy export from `stability._eigenvalues`). - [Structural sensitivity](core-concepts.md): `wavemaker` in `stability/sensitivity.py` combining direct/adjoint information for control and receptivity insight. - Factorization supports transpose solves needed by adjoint-related linear actions. ### Cross-workflow dependency summary ```mermaid flowchart LR CfgDom["ConfigurationManagement"] OrchDom["CaseOrchestration"] BFDom["BaseflowManagement"] MeshDom["MeshDiscretization"] OpDom["LNSEOperatorAssembly"] LADom["LinearAlgebra"] StabDom["StabilityAnalysis"] ResDom["ResolventAnalysis"] RtrDom["SolverRouting"] UtilDom["SharedUtilities"] CfgDom --> OrchDom OrchDom --> BFDom OrchDom --> OpDom OrchDom --> StabDom OrchDom --> ResDom OrchDom --> RtrDom OrchDom --> LADom BFDom --> MeshDom OpDom --> MeshDom OpDom --> BFDom OpDom --> LADom StabDom --> OpDom StabDom --> LADom ResDom --> OpDom ResDom --> LADom RtrDom --> OpDom RtrDom --> ResDom RtrDom --> LADom MeshDom --> UtilDom StabDom --> UtilDom ResDom --> UtilDom LADom --> UtilDom ``` ## Flow coordination and control ### Multi-module coordination mechanisms | Mechanism | How it works | |---|---| | **Staged pipeline ownership** | Driver owns stage order, logging, heartbeats, and payload; domains expose pure-ish services | | **Request/result objects** | Operator assembly uses structured requests (backend selection, sponge, finalization order) | | **Typed config seam** | `config/typed.py` mirrors `utils._solver_cfg` so call sites migrate without behavior change | | **Lazy public API** | Package `__init__` defers heavy imports; PETSc starts only on real solve paths (driver imports init early by design) | | **BC recipe orchestration** | `BCRule` / `BCRecipe` sequence declarative rules and return combined index sets for free DOFs | | **Router as campaign controller** | Map orchestration schedules β loops and frequency walks; LineWalker localizes reuse policy | | **Backend registry** | Central probe/normalize of factor backends before factorization commits resources | ### State management and synchronization **Logical case state progression:** ```mermaid stateDiagram-v2 [*] --> Unconfigured Unconfigured --> ConfigValid: schema_and_typed_merge ConfigValid --> BaseflowReady: load_and_preprocess BaseflowReady --> MeshReady: mesh_sbp_filters MeshReady --> OperatorReady: coeffs_L0_BCs OperatorReady --> FactorReady: factorize_shift FactorReady --> Solved: eigen_or_resolvent Solved --> Accepted: residual_acceptance_gates Accepted --> Emitted: payload_npz Emitted --> [*] ConfigValid --> Failed: ConfigValidationError FactorReady --> Failed: factor_backend_or_STRUMPACK_errors Solved --> Degraded: residual_gate_refactor_or_defer Degraded --> FactorReady: refactor_anchor Degraded --> Failed: walk_abort ``` **MPI/runtime synchronization:** - PETSc `COMM_WORLD` ranks share distributed matrices/vectors. - Rank 0 performs prefixed progress logs; stage complete lines are machine-parseable (`stage=... elapsed=... rss_peak_delta_bytes=...`). - Heartbeat thread/interval keeps long factors/EPS visible to operators without changing numerical state. - Factor objects encapsulate lifecycle: create → solve* N → refactorize (same nonzero pattern) → destroy. **Immutable/read-only defaults:** driver uses shared empty read-only NumPy arrays as hot-path defaults to avoid accidental mutation and allocation churn. ### Data passing and sharing | Data | Producer | Consumers | Form | |---|---|---|---| | Typed/merged config | config + yaml_io | all stages | dict / dataclasses | | BaseFlow | grid_fields + preprocess | coefficients, assembly, energy norm | dataclass, (ny,nx) arrays | | Mesh geometry/operators/runtime | mesh.py, sbp.py | assembly, BCs, filters | Mesh + view slices | | L0 | assembly.py | eigen, resolvent, router | SciPy CSC or PETSc Mat | | Free-DOF / BC index sets | BC recipe/dirichlet | reduced solves, validation | index arrays/masks | | Factorization handle | factorize.py | EPS shift-invert, H applies, walks | Factorization object | | Gains/modes/spectra | stability/resolvent | payload, map_record | NumPy + metadata dicts | | Cost samples | cost_metrics, router | certification, ops tuning | structured cost keys | | Provenance | provenance helpers | output naming, reproducibility | config MD5, git SHA, env | **Bridge policy:** SciPy ↔ PETSc conversion lives in `utils/petsc_core.py` to avoid circular imports between mesh and solver layers; complex-mode and distribution policy are encoded there. ### Execution control and scheduling 1. **Case-level scheduling**, driver selects eigen vs resolvent vs map stages from config; preflight advisories (config/baseflow/result) surface before/after heavy work. 2. **Eigen control**, shift-invert targets, backend from `solver.eigen` / env overrides; iterative eigen solver types listed explicitly in cfg helpers. 3. **Map scheduling**, - Outer loop: β - Inner: frequency line walk with anchor placement - Residual gate: continue / refactor / defer / abort - Optional stable sweep engine for smooth bands 4. **Resource control**, RAM autoselect, MUMPS ordering/BLR/precision options, optional GPU cuDSS path; cost break-even influences when reuse beats refactor. 5. **Concurrency model**, MPI process parallelism for distributed linear algebra; per-line anchor walks stay **sequential** to preserve factor reuse. Parallelism across β lines is an optimization opportunity when memory allows. 6. **Output control**, env overrides for config path, output NPZ, log prefix; atomic map writes reduce partial-result corruption. ## Exception handling and recovery ### Error detection and handling | Class | Detection | Typical response | |---|---|---| | **Config / schema** | Unknown keys, bad FD orders, invalid merged structure | Fail fast: `ConfigValidationError` before mesh/factor | | **Baseflow I/O** | Missing/invalid field names, non-2D floats | Raise at load; preflight baseflow advisories | | **BC recipe** | Invalid symbols, inconsistent borders/variables | Validation in recipe analyze/apply | | **Backend availability** | Requested MUMPS/LU/cuDSS not present | Probe in `backends.py`; error before factor commit | | **Factorization** | PETSc/MUMPS/STRUMPACK failures | Dedicated errors (e.g. STRUMPACK factorization failure path in eigen services); destroy partial state | | **Refactorize contract** | Nonzero pattern change | `NotImplementedError`/PETSc error if CSR pattern cannot be verified | | **Eigen/resolvent quality** | Residual norms, target distance, conjugate-pair gaps | Acceptance helpers; reorder/filter modes; advisories | | **Map walk quality** | Residual gate fail | Refactor anchor; defer walk; mark reason codes (OK/PENDING/FAILED) | | **Env parse** | Bad heartbeat seconds etc. | `ValueError` with clear message | ### Exception recovery mechanisms ```mermaid flowchart TD Detect["Detect failure or quality breach"] Class{"Failure class?"} FailFast["Fail fast free resources"] Refactor["Refactorize or new anchor"] Defer["Defer walk segment"] AcceptSoft["Soft accept with advisories"] RetryPol["Explicit policy retry only if configured"] EmitPartial["Emit certified partial map or error payload"] Abort["Abort case non-zero"] Detect --> Class Class -->|Config BC backend missing| FailFast Class -->|Residual gate map| Refactor Class -->|Walk budget exceeded| Defer Class -->|Mode quality marginal| AcceptSoft Class -->|Transient solver| RetryPol FailFast --> Abort Refactor --> EmitPartial Defer --> EmitPartial AcceptSoft --> EmitPartial RetryPol --> Detect EmitPartial --> Abort ``` **Recovery principles:** - **Prefer fail-fast on reproducibility gates** (config, missing backends) over running expensive wrong physics. - **Prefer refactor over silent inaccurate gains** on map residual gates. - **Iterative linear solvers are not an automatic downgrade** from direct factors; they require explicit solver policy (avoids hidden accuracy changes). - **Adjoint self-test / heal paths** exist in the factorization layer for iterative [KSP](solver-glossary.md) edge cases (bounded iterations), keeping adjoint consistency checks local to the LA layer. - **Resource cleanup:** factorization `destroy`, PETSc object lifecycle, and stage-scoped timing even on failure paths where implemented. ### Fault tolerance strategy design | Layer | Strategy | |---|---| | **Input** | Lenient-but-structured schema; preflight advisories without skipping hard validation | | **Numeric** | SBP dissipation policy mesh-aware; free-DOF reduction; residual validation on modes/gains | | **Linear algebra** | Direct default for reference accuracy; BLR/ordering knobs; optional GPU path; pattern-checked refactorize | | **Campaign** | Per-frequency reason codes; atomic savez; cost samples for post-mortem; certification metadata | | **Operations** | Heartbeats + stage RSS deltas for HPC job monitoring; rank-0 log prefix for aggregation | | **Reproducibility** | Config hash in outputs; git SHA provenance; env-isolated output paths for concurrent timing runs | ### Failure retry and degradation - **Retry:** Not a generic infinite retry loop. Map walks **locally** recover via refactor when residual gates fail. Explicit iterative policies may allow bounded KSP iterations. - **Degradation ladder (conceptual):** 1. Direct factor reuse (best cost/accuracy when residual OK) 2. Refactor at new anchor (restore accuracy, pay cost) 3. Stable-sweep / band-specific strategy 4. Backend change under RAM policy (e.g. memory-bound iterative, **opt-in**) 5. Defer/abort frequency or β line with recorded reason - **Partial success:** [Gain maps](core-concepts.md) can record completed (β, ω) cells with FAILED/PENDING markers rather than discarding the entire campaign when a single line breaks. - **Operator guidance:** On [STRUMPACK](https://portal.nersc.gov/project/sparse/strumpack/)/MUMPS failure, verify backend build, matrix [conditioning](https://en.wikipedia.org/wiki/Condition_number)/BC setup, and memory; do not "retry harder" without config change. ## Key process implementation ### Core algorithm processes #### LNSE temporal operator assembly **Intent:** Build large sparse \(L_0\) acting on \(q = [\rho, u, v, w, T]\) (size scale `5 * ngp`). **Pipeline characteristics (from assembly service):** 1. Prepare coefficients from BaseFlow (Cartesian vs cylindrical; axis conventions vs MATLAB port). 2. Attach FD/SBP differentiation and metrics. 3. Apply sponge and optional operator dissipation (\(L_0 \leftarrow L_0 - \varepsilon Q_{\mathrm{diss}}\) when enabled). 4. Enforce BC recipes ([Dirichlet](https://grokipedia.com/page/Dirichlet_boundary_condition) zero rows/cols + pin; [Neumann](https://grokipedia.com/page/Neumann_boundary_condition) as applicable). 5. Finalize backend: SciPy CSC and/or PETSc Mat with **fixed transformation order**. **Business rule:** Physics linearization is entirely driven by mesh-aligned mean flow; wrong preprocess (shock/sponge/viscosity) propagates directly into spectra and gains. #### Shift-invert eigenvalue process **Intent:** Extract eigenvalues of LNSE near targets via SLEPc EPS. **Technical spine:** 1. Form shifted operators as required by EPS SI strategy. 2. Factorize with selected backend (`resolve_factor_backend`, MUMPS options). 3. EPS iterations call repeated solves against the factor. 4. Post-process: residual norms, reordering, conjugate-pair checks, acceptance vs target. **Performance:** Factorization dominates; subsequent matvecs/solves amortize cost, same pattern as resolvent. #### Resolvent gain process **Intent:** Leading optimal gain and structures of \(H\). 1. Factor \((L_0 - (\gamma + i\omega) I)\). 2. Matvecs of \(H\) and \(H^{H}\) via B, solves, C (`ResolventOperator`). 3. Build `WeightedNormalOperator` shell for \(W\)-inner-product gain problem. 4. SLEPc on shell; validate residuals; return forcing/response pairs. **Coordination constraint:** SLEPc/MPI init ordering must be respected (shared with eigen path); package import stays light via type-only PETSc imports in places. #### Anchor-preconditioned frequency walk **Intent:** Minimize factor count across dense ω grids. ```mermaid flowchart TD WalkInit["LineWalker start band"] PayAnchor["Factor true operator at anchor omega"] MeasureCost["Record factor time memory compression"] NextOm["Next frequency"] PCApply["Use anchor factor as right PC"] SolveOm["Solve shifted system gain EPS"] Gate{"Residual below threshold?"} Keep["Accept gain continue"] NewFac["New real factor become anchor"] Done{"Line complete?"} WalkEnd["Return samples reasons"] WalkInit --> PayAnchor PayAnchor --> MeasureCost MeasureCost --> NextOm NextOm --> PCApply PCApply --> SolveOm SolveOm --> Gate Gate -->|Yes| Keep Gate -->|No| NewFac NewFac --> MeasureCost Keep --> Done Done -->|No| NextOm Done -->|Yes| WalkEnd ``` ### Data processing pipelines #### Mean-flow to operator pipeline ``` CFD npy/npz → field validation (names, shapes) → BaseFlow dataclass → shock cleanup / sponge / Sutherland (or configured) viscosity → Mesh geometry (+ deformation, masking, poles) → SBP D, norms, filters, optional Q_diss → physical coefficient blocks → sparse assembly + sponge + dissipation → BC enforcement + free-DOF sets → SciPy and/or PETSc L0 ``` #### Result and metrics pipeline ``` Solver raw outputs → residual/acceptance filters → payload builder / map_record → cost_metrics (wall stages, RSS, GPU mem, ranks) → provenance (config hash, git SHA, PETSc options tags) → NPZ / atomic savez emission ``` Driver stage logging provides an ops-facing timeline parallel to scientific outputs. ### Business rule execution | Rule | Enforcement point | |---|---| | Reproducible cases are YAML-driven | config merge + schema before heavy stages | | State ordering is ρ,u,v,w,T | assembly and public API docs | | Direct sparse solve is default accuracy path | factorize policy; tests/examples | | Iterative LA only when explicitly chosen | solver-type allowlists | | Map dissipation needs dual enable | router docs + assembly strength read from top-level `dissipation_epsilon` | | Spatial config may inherit eigen policy | spatial_eig fallback | | Library import must stay side-effect light | lazy `_LAZY_EXPORTS` in package `__init__` | | Rank-0 logging only for human progress | `_log` guards | | Gains/modes require residual accountability | analysis validation + acceptance modules | ### Technical implementation details #### Runtime driver (`driver/run.py`) - Initializes PETSc from argv at import of the driver module (case path intentionally heavy). - Resolves config/output paths via environment with documented fallbacks. - Stages work with perf counters and peak RSS deltas for HPC observability. - Wires: `load_case_fields`, `Mesh`, `build_lnse_system_result`, BC analysis, `temporal_eigenvalues` / `_solve_near_target`, `resolvent_analysis`, router helpers, `build_runtime_payload`. - Stability acceptance and conjugate-pair checks gate scientific quality before trusting payloads. #### Factorization layer (`linalg/factorize.py`) - Factorize-once, solve-many for SI LNSE workflows. - Default **KSPPREONLY + PCLU + MUMPS**. - Supports solve, transpose solve, block solve, refactorize, destroy. - Optional cuDSS factor path as PC shell. - Lazy PETSc accessor pattern so non-driver imports need not initialize PETSc until use. #### Router (`solvers/router`) - `resolvent_map`: turnkey β/frequency orchestration; largest hub in package. - `map_walk`: measured anchors + residual gate. - `stable_sweep`: anchor-preconditioned stable-regime engine. - Supporting: anchors, recycling, break_even, ram_autoselect, certification, band_marking, map_record. #### Mesh/SBP (`mesh`) - `Mesh` with geometry / operators / runtime-state views for explicit consumer contracts. - `sbp.py`: Strand / Mattsson–Nordström style operators; pole matrices for axisymmetric grids. #### External runtime contract | External system | Role in workflows | |---|---| | PETSc | Distributed Mat/Vec, KSP/PC | | SLEPc | EPS for eigen and gain problems | | MUMPS | Default sparse direct factor | | MPI | Process parallelism | | cuDSS | Optional GPU direct factor | | NumPy/SciPy | I/O, CSC/CSR assembly, local kernels | | YAML + npy/npz | Sole primary I/O surface for cases and baseflow | ### Performance optimization processes | Strategy | Where applied | Effect | |---|---|---| | Factorize-once solve-many | Eigen SI, resolvent H, map PC reuse | Amortizes dominant cost | | Anchor + residual gate | Frequency lines | Fewer full LU factors at controlled accuracy | | Stable sweeps | Smooth/stable bands | Cheap ω coverage with preconditioners | | Backend/RAM autoselect | Large cases | Avoids OOM; chooses MUMPS/BLR/cuDSS/iterative policy | | SciPy assemble → PETSc solve | Operator path | Fast local build, scalable solve | | Lazy imports | Library users | Fast startup; PETSc only when needed | | Stage RSS/wall metrics | Driver + router | Data-driven anchor and backend tuning | | Optional dissipation off on rectilinear meshes | Map prep | Avoids nnz/time penalty with negligible gain change | | MPI distribution | Global Mat/Vec | Memory scale-out for large ngp | **Optimization opportunities (architecture-aligned):** 1. Deeper reuse of one prepared L0 across mixed temporal + resolvent stages in a single case. 2. Residual-gate thresholds tuned per backend accuracy profile (MUMPS BLR vs cuDSS). 3. Parallel β lines when memory permits; keep per-line walks sequential for factor reuse. 4. Richer structural-sensitivity packaging on the standard payload path. 5. Feed multi-node cost samples more aggressively into anchor placement. ### Operability guidance **Development** - Extend physics in coefficients/assembly; keep BC recipes declarative. - New analysis products should consume finalized L0 + factorization APIs rather than bypassing them. - Preserve lazy export surface in `__init__.py` for embedders. **Operations / HPC** - Monitor `[prefix] stage complete: ...` and heartbeats for stuck factors. - On map jobs, track `cost_wall_factor_s` vs `cost_wall_solve_s` and residual-driven refactor rate. - Set `LINSTABPY_CASE_CONFIG` and dedicated `LINSTABPY_CASE_OUTPUT_NPZ` per concurrent run. - Verify dissipation flags in pairs on curvilinear meshes; confirm undamped intent on rectilinear. **Troubleshooting quick map** | Symptom | Likely stage | First checks | |---|---|---| | Immediate abort | Config/schema | YAML keys, FD orders | | Bad spectra/gains | Baseflow/BC/assembly | Field alignment, BC recipe, sponge | | Factor crash / OOM | Linear algebra | Backend probe, ranks, BLR, iterative opt-in | | Noisy map ridges | Router/mesh dissipation | metric cross-coupling, epsilon dual enable | | Slow dense ω grid | Map walk | Anchor spacing, residual threshold, stable sweep | ## Document summary LinStabPy core workflows form a **single staged LNSE pipeline** with multiple scientific products. **Baseflow ingest and operator prep** produce a BC-enforced \(L_0\); **temporal/spatial/adjoint eigenanalysis** and **resolvent analysis** share that operator and a **factorize-once** PETSc/MUMPS spine; **β–ω gain maps** add cost-aware routing with anchors, preconditioner reuse, and residual gates. Coordination is hub-and-spoke through the case driver and router, with configuration as the reproducibility gate and residual/acceptance checks as the scientific quality gate. Exception handling favors fail-fast validation, explicit accuracy policies, and campaign-local refactor/defer rather than silent degradation, supporting reliable MPI-scale stability and resolvent studies for research and aerospace engineering users. --- ## Module guide Source: https://linstabpy.vercel.app/modules/ # Module guide Short reference per package, then longer notes. Policy lives in the owning module; `docs/architecture.md` is the short map. ## Quick reference | Package | Page | |---|---| | `driver` | [Runtime driver](driver.md) | | `operators` | [LNSE assembly](operators.md) | | `boundary_conditions` | [BC recipes](boundary-conditions.md) | | `mesh` | [Structured mesh](mesh.md) | | `baseflow` | [Packed fields](baseflow.md) | | `config` | [YAML and Config](config.md) | | `stability` | [Eigenvalues](stability.md) | | `resolvent` | [Resolvent](resolvent.md) | | `solvers` | [Solvers](solvers.md) | | `linalg` | [Sparse factors](linalg.md) | ## Package notes | Topic | Page | |---|---| | Case pipeline | [Case Orchestration](case-orchestration.md) | | Configuration | [Configuration Management](configuration-management.md) | | Baseflow | [Baseflow Management](baseflow-management.md) | | Mesh | [Mesh Discretization](mesh-discretization.md) | | LNSE assembly | [LNSE Operator Assembly](lnse-operator-assembly.md) | | Stability | [Stability Analysis](stability-analysis.md) | | Resolvent | [Resolvent Analysis](resolvent-analysis.md) | | Solver routing | [Solver Routing and Sweep Optimization](solver-routing.md) | | Linear algebra | [Linear Algebra Infrastructure](linear-algebra.md) | | Utilities | [Shared Scientific Utilities](shared-utilities.md) | --- ## Solver glossary Source: https://linstabpy.vercel.app/solver-glossary/ # PETSc and SLEPc glossary Definitions of every solver term used in this wiki, linked to the official docs. Assumes linear algebra ([Krylov](https://grokipedia.com/page/Krylov_subspace) subspaces, [LU](https://grokipedia.com/page/LU_decomposition), [eigenvalues](https://grokipedia.com/page/Eigenvalues_and_eigenvectors)) No prior PETSc or SLEPc experience needed. Readers who already use PETSc can skip to [The exact chains this repo runs](#the-exact-chains-this-repo-runs). **The two libraries.** [PETSc](https://petsc.org/release/manual/ksp/) (Portable, Extensible Toolkit for Scientific Computation) owns distributed vectors and matrices plus linear solvers. [SLEPc](https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetProblemType/) (Scalable Library for Eigenvalue Problem Computations) sits on top of PETSc and owns eigensolvers. LinStabPy drives both through their Python bindings, [petsc4py](https://petsc4py.readthedocs.io/en/stable/) and [slepc4py](https://slepc4py.readthedocs.io/en/stable/), over [mpi4py](https://mpi4py.readthedocs.io/) for [MPI](https://grokipedia.com/page/Message_Passing_Interface). The install guide explains how to build this stack. This page explains what the objects mean. Physics and math vocabulary (growth rate, gain, Chu norm, wavemaker) lives in the [core concepts](core-concepts.md) companion. Repo conventions used below: $\omega = \mathrm{eig}(i L_0)$, unstable exactly when $\mathrm{Im}(\omega) > 0$; [resolvent](https://grokipedia.com/page/Resolvent_formalism) $H = C F (L_0 - (\gamma + i\omega) I)^{-1} F B$. ## Mental model Four object pairs cover nearly everything: | Pair | Library | Job in one line | |---|---|---| | [`Mat`](https://petsc.org/release/manualpages/Mat/Mat/) and [`Vec`](https://petsc.org/release/manualpages/Vec/Vec/) | PETSc | A sparse matrix and a vector, distributed across MPI ranks. | | [`KSP`](https://petsc.org/release/manual/ksp/) and [`PC`](https://petsc.org/release/manual/ksp/) | PETSc | A linear solver (`KSP`) split from its preconditioner (`PC`), so any iteration can pair with any preconditioner. | | [`EPS`](https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetProblemType/) and [`ST`](https://slepc.upv.es/documentation/current/docs/manualpages/ST/STSetShift/) | SLEPc | An eigensolver (`EPS`) split from its spectral transformation (`ST`), so any eigensolver can pair with shift-invert. | | Shell matrices | PETSc | A `Mat` you never assemble. You provide the multiply function and PETSc calls it. | A typical LinStabPy solve wires all four: a distributed `Mat` holds the shifted operator, a `KSP` with `PC` factors it once, an `EPS` with `ST` runs shift-invert Krylov-Schur calling that factor on every iteration, and the resolvent path adds a shell matrix whose multiply applies $H$ through the same factor. ## PETSc basics **PETSc** is a C library (with Python bindings) for vectors, matrices, and solvers on distributed memory. Start with the [KSP manual chapter](https://petsc.org/release/manual/ksp/). It covers linear solves end to end and is the best single document behind this page. **Communicators.** Every PETSc object lives on an MPI communicator, almost always `COMM_WORLD` here. Ranks own disjoint row ranges of each matrix and vector. Only the owning rank writes its range. When the docs say "rank 0 owns logging and IO decisions," that is convention in our driver, not a PETSc requirement. **[`Mat`](https://petsc.org/release/manualpages/Mat/Mat/) and [`Vec`](https://petsc.org/release/manualpages/Vec/Vec/).** The distributed matrix and vector types. The default [sparse](https://grokipedia.com/page/Sparse_matrix) CPU format is AIJ (compressed sparse row across ranks). LinStabPy assembles locally in SciPy CSC or CSR, then converts to a distributed `Mat` in `utils/petsc_core.py`. Sizes here are $5$ unknowns per grid point before boundary elimination, so a $300 \times 100$ grid gives a $150{,}000 \times 150{,}000$ complex sparse matrix. **SciPy CSC and CSR.** The two compressed formats from [`scipy.sparse`](https://docs.scipy.org/doc/scipy/reference/sparse.html) (background: [sparse matrix](https://grokipedia.com/page/Sparse_matrix)). CSR stores row by row: values plus column indices plus row-start pointers. CSC stores the transpose pattern: values plus row indices plus column-start pointers. Same entries, different access direction. Assembly here builds small blocks in CSR and finalizes the full $L_0$ as complex CSC. Both convert to each other freely, and every solver entry point accepts either; the bridge normalizes to CSR on the way into PETSc. **PETSc AIJ.** PETSc's sparse format: a distributed CSR. Each rank owns a contiguous row range and stores its rows as `(indptr, indices, data)` triples, created via `Mat().createAIJ(size, csr, comm)`. Two layouts exist. *Local* AIJ on `COMM_SELF` holds the whole matrix on one rank (single-rank style). *Distributed* AIJ on `COMM_WORLD` splits rows across ranks by ownership ranges. Same object type, different communicator and row distribution. **The bridge (`utils/petsc_core.py`).** One module owns every conversion so call sites never hand-roll them: | Function | Job | |---|---| | `scipy_to_petsc` | SciPy matrix to PETSc AIJ `Mat`. Layout follows `petsc_matrix_mode` (below). | | `scipy_to_petsc_distributed`, `scipy_to_petsc_distributed_from_root` | Ownership-aware distributed AIJ, scattering rows to owner ranks. | | `petsc_to_scipy` | PETSc `Mat` back to SciPy (`csr` or `csc`). Used for comparison, parity, and SciPy-only paths. | | `create_distributed_aij`, `insert_local_csr_rows`, `assemble_matrix` | Lower-level distributed construction used inside assembly and conversion. | | `ownership_range`, `petsc_owned_rows` | Which global rows this rank owns, and their CSR data. | | `require_complex` | Asserts complex-scalar PETSc and installs one-time MUMPS defaults. Fails fast otherwise. | | `compare_sparse_matrices`, `compare_petsc_against_scipy`, `petsc_matvec_residual` | Parity helpers proving the two paths hold the same operator. | **`assembly_backend`: where $L_0$ gets built.** YAML or `LINSTABPY_ASSEMBLY_BACKEND`, default `scipy`. `scipy` assembles CSC locally and converts on demand: the cross-machine baseline, and mandatory for adjoint workflows (`adjoint_eigenvalues` rejects PETSc matrices today). `petsc` assembles AIJ directly, but with two documented limits: direct PETSc assembly cannot apply `dissipation_epsilon`, and distributed direct assembly covers model `2D` only. When in doubt, SciPy-assemble then convert. **`petsc_matrix_mode`: local or distributed AIJ.** YAML or `LINSTABPY_PETSC_MATRIX_MODE`, default `auto`. `local` and `serial` force single-rank AIJ. `distributed`, `mpi`, and `parallel` force ownership-split AIJ. `auto` picks distributed exactly when the MPI world has more than one rank (and local when PETSc is absent). Anything else raises `ValueError` naming the valid words. **Rule of thumb.** Small cases, debugging, adjoint work, and MATLAB parity: `scipy` backend with local matrices. Production MPI runs: same SciPy assembly, converted to distributed AIJ for the factorize and eigensolve stages. The operator values are identical either way; only the layout changes. **Complex scalars.** PETSc compiles for either real or complex floating point, chosen at build time. This repo requires complex builds: eigenvalues, shifts, and modes are complex throughout. `install.sh --check` asserts this. A real-scalar PETSc imports fine and then fails real work, which is why the check exists. **Init order.** Import `mpi4py` before `petsc4py` so MPI initializes with thread level `MULTIPLE`, which some factor backends need. Plain `import linstabpy` initializes neither. The driver and solver modules initialize them when a solve stage starts. Solver modules that need SLEPc call the shared init helper first. **Options database.** PETSc and SLEPc read command-line style options at runtime (`-ksp_type`, `-eps_monitor`, and hundreds more). Any example can receive them after `--` through the runner, for example `-- -eps_monitor` to watch eigensolver convergence. YAML and env policy in this repo set the important ones explicitly, but the options database is always available as an override and debugging channel. ## KSP: linear solvers and the PREONLY direct path **[`KSP`](https://petsc.org/release/manual/ksp/)** (Krylov Subspace Solver) is PETSc's linear solver object: given $Ax = b$, it produces $x$. The key design point is the split: `KSP` chooses the *iteration* (or none), [`PC`](https://petsc.org/release/manual/ksp/) chooses the *[preconditioner](https://grokipedia.com/page/Preconditioner)*. You can pair [GMRES](https://grokipedia.com/page/Generalized_minimal_residual_method) with LU, or no iteration with LU, without touching calling code. SLEPc shift-invert and our resolvent applies both consume `KSP` objects through this interface. **[`KSPPREONLY`](https://petsc.org/release/manualpages/KSP/KSPPREONLY/)** means: no Krylov iteration at all, just apply the preconditioner once. When the preconditioner is an exact LU factorization, that single apply *is* the direct solve. This is why our default shifted solve is described as "`KSPPREONLY` with `PCLU` through [MUMPS](https://en.wikipedia.org/wiki/MUMPS_(software))": mathematically it is sparse direct LU, reached through the `KSP` interface so every consumer (eigen shift-invert, resolvent applies, transpose solves) shares one code path. Stated in `linalg/factorize.py`: the default path is `KSPPREONLY` + `PCLU` + MUMPS. **[`KSPGMRES`](https://petsc.org/release/manualpages/KSP/KSPGMRES/)** is restarted GMRES, the standard Krylov method for nonsymmetric systems. GMRES builds an orthogonal basis of the Krylov subspace and minimizes the residual over it, restarting every *restart* steps to bound memory (default 50 in our router). This repo uses GMRES in exactly one place: the $\beta$ and $\omega$ map router reuses an anchor factorization as a right preconditioner for nearby frequencies, solving each nearby system with right-preconditioned GMRES (restart 50, max 200 iterations) until a residual gate forces a fresh factorization. See the Solver Routing page. **[`KSPFGMRES`](https://petsc.org/release/manualpages/KSP/KSPFGMRES/)** is flexible GMRES, which tolerates a preconditioner that changes between iterations (for example an inner iterative solve). Ours uses it on two explicit opt-in paths: `fgmres_lu` (FGMRES around an LU preconditioner) and `strumpack_fgmres` (FGMRES around a [STRUMPACK](https://portal.nersc.gov/project/sparse/strumpack/) approximate factor), plus the [cuDSS](https://docs.nvidia.com/cuda/cudss/) GPU path where the GPU factor preconditions an FGMRES outer iteration. Iterative paths are never silent fallback here. They require explicit solver policy. **[`KSPSetTolerances`](https://petsc.org/release/manualpages/KSP/KSPSetTolerances/)** sets relative tolerance, absolute tolerance, divergence tolerance, and maximum iterations for the Krylov iteration. On the default `KSPPREONLY` path these knobs do nothing, because there is no iteration. Accuracy there comes from the factorization quality and from EPS-side settings instead. The `ksp_rtol`, `ksp_atol`, and `ksp_max_it` knobs matter only on iterative paths, which is a common source of confusion worth stating twice. ## PC: preconditioners, LU, and MUMPS **`PC`** (preconditioner) is the other half of the `KSP` split. A preconditioner approximates $A^{-1}$ cheaply so the outer iteration converges fast, or in our default case *is* $A^{-1}$ computed exactly by factorization. **[`PCLU`](https://petsc.org/release/manualpages/PC/PCLU/)** means the preconditioner is an LU factorization: on setup it factors $A = LU$ (with pivoting), on apply it runs triangular solves. With `KSPPREONLY` in front, setup plus one apply equals a direct solve. In parallel the factorization itself is distributed (see MUMPS below). **[`PCFactorSetMatSolverType`](https://petsc.org/release/manualpages/Mat/MATSOLVERMUMPS/)** selects *which* sparse direct package performs the LU inside `PCLU`. This one call is our entire backend seam. Resolution order in this repo: YAML `factor_backend`, then per-solver env override, then `LINSTABPY_FACTOR_BACKEND`, then the call-site default. An unavailable backend fails loudly during setup. It never falls back silently. **[`MATSOLVERMUMPS`](https://petsc.org/release/manualpages/Mat/MATSOLVERMUMPS/) ([MUMPS](https://mumps-solver.org/))** is the default backend: a [MUltifrontal](https://en.wikipedia.org/wiki/Multifrontal_method) Massively Parallel sparse direct Solver that factors distributed matrices across MPI ranks. "Multifrontal" names the algorithm family (dense frontal matrices eliminated in a tree order). Practical knobs exposed here: ordering selection mapped to MUMPS `ICNTL(7)` (`amd`, `pord`, `metis`, `qamd`, `auto`), optional block low-rank compression (`LINSTABPY_MUMPS_BLR` with tolerance), and precision selection. MUMPS uses dynamic pivoting, which is why the temporal eigensolver retries a failed STRUMPACK factor (static pivoting, `ZERO_PIVOT`) once with MUMPS. **Other backends, one line each.** STRUMPACK: structured and rank-compressed direct and approximate solver, experimental here, used through FGMRES. [SuperLU_DIST](https://portal.nersc.gov/project/sparse/superlu/): distributed-memory sparse LU, experimental here, GPU offload unvalidated. cuDSS: NVIDIA GPU sparse direct LU, validated here *as a preconditioner inside FGMRES* (one rank, one GPU, with a second factorization of the conjugate transpose for [adjoint](https://grokipedia.com/page/Adjoint) applies), matching MUMPS to 1e-12 on tested cases, not as a PETSc-native factor type. The backend audit lives in `docs/COMPUTE_BACKENDS.md`. **[`PCBJACOBI`](https://petsc.org/release/manualpages/PC/PCBJACOBI/) and [`PCASM`](https://petsc.org/release/manualpages/PC/PCASM/)** are [domain-decomposition](https://grokipedia.com/page/Domain_decomposition_methods) preconditioners: block Jacobi (each rank factors its own diagonal block, no overlap) and additive Schwarz (overlapping blocks). They appear as the baseline direction for the memory-bound iterative eigensolver path. One API trap the code documents: sub-preconditioner options must be applied *after* `EPS.setUp()`, or they attach to nothing. **Transpose solves.** Adjoint work and $H^H$ applies need $A^H x = b$, not $Ax = b$. The factorization object exposes `solve_transpose` alongside `solve`, reusing the same LU factors. No second factorization is needed for the adjoint half of resolvent or stability analysis. ## MatShell: matrices you never assemble **[`MatCreateShell`](https://petsc.org/release/manualpages/Mat/MatCreateShell/)** creates a PETSc `Mat` with no stored entries. You register callback functions (at minimum the multiply `y = Ax`, optionally transpose and others) and PETSc treats the result as a matrix: it can sit inside `KSP`, `EPS`, and `ST` exactly like an assembled `AIJ` matrix. "Shell" contrasts with "assembled." The multiply runs your code, typically a few operator applies and triangular solves. **Why the resolvent needs it.** The resolvent $H = C F (L_0 - (\gamma + i\omega) I)^{-1} F B$ is dense even though $L_0$ is sparse. Forming it explicitly at $150{,}000$ unknowns is out of the question. Instead `WeightedNormalOperator` is a shell whose multiply chains existing cheap pieces: apply $H$ (two filters plus one factored solve plus $B$ and $C$ maps), weight by $W$, apply $H^H$ (adjoint filters plus one transpose solve), and for the [Chu](core-concepts.md) forcing norm multiply by $W^{-1}$. SLEPc then runs an eigensolver that only ever calls this multiply. [Gains](core-concepts.md) come out as $\sigma = \sqrt{\mathrm{Re}(\lambda)}$ without $H$ ever existing as entries. **Ownership detail.** On multi-rank communicators the multiply callback writes only its rank's local row range into the result vector. Getting this wrong silently corrupts outer iterations, which is why the shell code is explicit about ranges. ## SLEPc: EPS eigensolvers and ST spectral transformations **SLEPc** extends PETSc to eigenvalue problems. It operates on PETSc `Mat` and `Vec` objects, so everything in sections 2 to 5 carries over. Its Python binding is [slepc4py](https://slepc4py.readthedocs.io/en/stable/). **[`EPS`](https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetProblemType/)** (Eigenvalue Problem Solver) is the eigensolver object, the analogue of `KSP` for $Ax = \lambda Bx$. You set the problem type, the method, how many pairs you want, which ones, the tolerance, and the spectral transformation, then call solve and extract pairs one by one. **Problem types ([`EPSSetProblemType`](https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetProblemType/)).** This repo uses two. `NHEP`: standard ($B = I$) non-Hermitian problem, used for temporal stability ($\omega = \mathrm{eig}(iL_0)$) and for the weighted normal resolvent problem. `GNHEP`: generalized non-Hermitian $Ax = \lambda Bx$, used for spatial stability where the $B$ matrix carries boundary rows and algebraic constraints (and is often singular, which matters below). Hermitian types exist and are not used here because LNSE operators are [nonnormal](https://en.wikipedia.org/wiki/Defective_matrix). **[`EPSKRYLOVSCHUR`](https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSKRYLOVSCHUR/)** is the restarted Krylov-Schur eigensolver, our method on every path. Think restarted [Arnoldi](https://grokipedia.com/page/Arnoldi_iteration) with a cleaner restart: build a Krylov subspace, project the problem, keep the wanted Ritz information in a Krylov-Schur decomposition, expand again. After $m$ Arnoldi steps the factorization reads $AV = VH + \beta v_{m+1} e_m^T$ with $V$ orthonormal and $H$ upper Hessenberg. Each eigenpair $(\lambda_i, y_i)$ of the small $H$ lifts to a **Ritz pair** $(\lambda_i, Vy_i)$ of $A$, with residual exactly $|\beta|\,|e_m^T y_i|$, so convergence is tested with no extra matvecs. Krylov-Schur (Stewart, 2001) keeps $H$ in Schur form instead, which makes restarting clean: truncate to the wanted Schur vectors and expand again, bounding memory at subspace size rather than iteration count. Like all Krylov eigensolvers it finds exterior eigenvalues of the operator it sees, which is exactly why shift-invert (next object) precedes it for interior targets. Full development in Frantz, Loiseau and Robinet (2023), *Krylov methods for large-scale dynamical systems*, sections on Krylov subspaces and large-scale eigensolvers. **[`EPSSetDimensions`](https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetDimensions/)** sets `nev` (eigenpairs wanted) and `ncv` (Krylov subspace size, the accuracy and cost knob). `ncv` must exceed `nev`, and larger `ncv` helps crowded or nonnormal neighborhoods where small subspaces stall. Repo defaults: temporal uses $\max(2\cdot\mathrm{nev}+3, 30)$ clipped to matrix size; spatial uses $\min(\max(4\cdot\mathrm{nev}, 30), \mathrm{ndof}-1)$; resolvent uses $\min(\max(2\cdot\mathrm{nev}+3, 20), \mathrm{ndof}-1)$. A user-passed `ncv` must be larger than `nev` or setup raises. **[`EPSSetWhichEigenpairs`](https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetWhichEigenpairs/)** picks *which* eigenvalues to want. Temporal and spatial paths use `TARGET_MAGNITUDE`: closest to the shift target in magnitude, combined with shift-invert. The resolvent path uses `LARGEST_MAGNITUDE`: the normal operator's largest eigenvalues are the squared leading gains by construction, so no shift is needed there. **Krylov convergence practice.** Two results from the review govern expectations. First, isolated eigenvalues converge geometrically fast while branches crawl: the error bound scales with $\kappa(V)\,(\rho/|\mu_i - c|)^k$ for an eigenvalue $\mu_i$ outside a disk of radius $\rho$, and dense branches are the worst case on both factors (poor isolation times large eigenvector conditioning from nonnormality). Crowded or branch-heavy targets are exactly where raising `ncv` pays. Second, Arnoldi favors eigenvectors with large components in the starting vector, so genuinely random starts are the default unless prior knowledge says otherwise; a start nearly orthogonal to the wanted mode stalls. Time-stepper codes add heuristics this repo does not need (matrix-forming factors replace sampling periods), but the diagnosis transfers: slow convergence near a branch means widen the subspace or move the shift, not just tighten the tolerance. See Frantz, Loiseau and Robinet (2023), observations on eigenvalue convergence. **[`EPSSetTolerances`](https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetTolerances/)** sets the eigenpair tolerance and restart limit. Three finer points the code exposes: the *convergence test* (`abs`, `rel`, or `norm`) decides how the residual is measured during the solve; the *reported residual* need not equal the tested one; and `true_residual=true` forces the expensive untransformed residual on shift-invert runs, because the cheap transformed estimate can look converged while the true pair is poor. Temporal and resolvent tolerances come from YAML (`eps_tol`, `eps_max_it`); spatial currently uses fixed `tol=1e-10`, `max_it=max(500, 10\cdot\mathrm{ndof})`. Monitor output (`-eps_monitor`) prints per-restart convergence and is the first thing to ask for when a spectrum looks wrong. **[`ST`](https://slepc.upv.es/documentation/current/docs/manualpages/ST/STSetShift/) and [`STSINVERT`](https://slepc.upv.es/documentation/current/docs/manualpages/ST/STSINVERT/)** (spectral transformation, shift-invert). Krylov methods find exterior eigenvalues, but stability targets sit in the spectrum interior. Shift-invert replaces the operator with $(A - \sigma B)^{-1} B$, whose *largest* eigenvalues are exactly the $\lambda$ closest to the shift $\sigma$, with eigenvectors unchanged. Each `ST` apply is one linear solve with the shifted matrix, which is why the `KSPPREONLY` + `PCLU` + MUMPS factor from sections 3 and 4 sits underneath every shift-invert eigensolve. [`STSetShift`](https://slepc.upv.es/documentation/current/docs/manualpages/ST/STSetShift/) sets $\sigma$. This repo sets the shift explicitly on `ST` rather than relying only on `EPS.setTarget`, so the generalized shift center is never implicit. Temporal targets live in $\omega$ space ($\sigma$ near the wanted $\omega$); the solver internally handles the $iL_0$ scaling. ## The exact chains this repo runs One scoping note: this repo is matrix-forming (assembled $L_0$ plus direct factors), not matrix-free time-stepping, so the review's Newton-Krylov, monodromy, and sampling-period machinery does not apply here. What transfers fully is the Krylov eigensolver core below: Arnoldi factorization, Ritz convergence, and shift-invert. | Analysis | Operator seen by EPS | EPS setup | ST and linear backend | |---|---|---|---| | **Temporal EVP** | $iL_0$, target $\sigma$ in $\omega$ space | `NHEP`, Krylov-Schur, `TARGET_MAGNITUDE`, YAML tol and max it | `SINVERT` at $\sigma$, factorize-once `KSPPREONLY` + `PCLU` + MUMPS | | **Spatial EVP** | $(A, R)$ block pair | `GNHEP`, Krylov-Schur, `TARGET_MAGNITUDE`, `tol=1e-10` | `SINVERT`, same direct backend. Tiny or near-full problems (`ndof ≤ 2` or `nev ≥ ndof-1`) switch to dense QZ (`scipy.linalg.eig`) because $R$ is singular from BC rows and QZ handles that natively. | | **Resolvent gain** | Shell: $W^{-1}H^HWH$ (`chu`) or $H^HWH$ (`l2`) | `NHEP`, Krylov-Schur, `LARGEST_MAGNITUDE`, YAML tol and max it | No `ST`. Each shell multiply applies $H$ and $H^H$ through one `KSPPREONLY` + `PCLU` factor of $L_\sigma$. | | **Router reuse** | Nearby $K(\omega)$ | Same gain EPS per point | Anchor factor as *right preconditioner* inside GMRES (restart 50, max 200). Residual gate decides reuse, refactor, or defer. | | **Adjoint EVP** | $W^{-1}L_0^H W$ (SciPy path) | Same temporal EPS shape near mapped target | Transpose solves through the existing factors. | | **Iterative opt-ins** | Same operators | Same EPS shape | `fgmres_lu` or `strumpack_fgmres`, explicit policy only, with `true_residual` defaulted on because transformed residuals can lie. | Two boundaries worth memorizing. A direct factor belongs to exactly one shift $(\omega, \gamma)$. Reuse across shifts happens only as router preconditioning under a residual gate, never as an exact solve. And iterative linear solvers are never automatic downgrade from direct factors. They require explicit solver policy so accuracy changes stay visible. ## Knobs: YAML, env, and options database | Knob | Where | Effect | |---|---|---| | `solver.eigen` / `solver.spatial` / `solver.resolvent` YAML blocks | Case file | Targets, `nev`, `ncv`, `eps_tol`, `eps_max_it`, convergence test, true residual flag, monitor flag per analysis. Spatial falls back to the eigen block when absent. | | `factor_backend` | YAML solver subsection, then `LINSTABPY_FACTOR_BACKEND`, then call-site default | Which package backs `PCLU`: `mumps` (default), `strumpack`, `superlu_dist`, `cudss` route. Per-solver env overrides exist for eigen and resolvent paths. | | `LINSTABPY_MUMPS_ORDERING`, `LINSTABPY_MUMPS_BLR`, `LINSTABPY_MUMPS_BLR_TOL`, `LINSTABPY_FACTOR_PRECISION` | Env | MUMPS ordering, low-rank compression, factor precision on paths that consume them. | | `ncv`, `eps_tol`, `eps_max_it`, `conv_test`, `true_residual`, `monitor` | YAML `eps` subsection | Subspace size, tolerance, restart cap, how convergence is measured, whether to compute true residuals, whether to print monitor lines. | | `-eps_monitor`, `-ksp_monitor`, `-eps_view`, … | Options database (after `--`) | Runtime inspection without editing files. `-eps_monitor` is the single most useful flag when a spectrum misbehaves. | | `ksp_rtol`, `ksp_atol`, `ksp_max_it`, `pc_type` | YAML, iterative paths only | Meaningful only when a real Krylov iteration runs. Inert on `KSPPREONLY`. | ## Choosing a configuration Start here when the knobs in [Knobs](#knobs-yaml-env-and-options-database) feel like too many degrees of freedom. Every number below is measured on the Savarino resolvent case (450k unknowns, 23.4M nonzeros, dual Xeon Silver 4210R) unless stated otherwise. ```mermaid flowchart TD Q1{"Does a direct factor fit in memory?"} Q1 -->|Yes| D1["Direct: KSPPREONLY + PCLU + MUMPS"] Q1 -->|No| I1["Iterative, explicit opt-in: fgmres_lu or strumpack_fgmres with true_residual on"] D1 --> Q2{"One point or a map?"} Q2 -->|One point| P1["Single factorize, single EPS. Done."] Q2 -->|Map over frequencies| Q3{"Lines or farm?"} Q3 -->|Few points, many nodes| F1["Farm: independent jobs per point (83.6 min for 200 points, 5 concurrent)"] Q3 -->|Dense lines| L1["Line driver: assemble once, anchor factors, in-place numeric refactorization (59.7 min for the same 200)"] D1 --> Q4{"Which mesh?"} Q4 -->|Rectilinear, leak at round-off| E0["epsilon = 0. Check insensitivity under refinement."] Q4 -->|Bent, leak order one| E1["Dissipation on. Tune epsilon on this mesh until the gain or eigenvalue stops moving."] ``` Why the defaults are the defaults, with the measurements behind them: | Decision | Evidence | |---|---| | Direct first | One point costs 60 s factorize from scratch (54 s numeric refactor reusing symbolic analysis), 0.44 s per forward solve, 0.43 s per adjoint solve, 5.8 s for the gain eigensolve. The factor is four fifths of a point; everything else is noise next to it. | | Small Krylov spaces | Cutting `ncv` from 20 to 5 cut solves per point from 40 to 10 and moved the gain by less than $10^{-10}$. Subspace size sets cost, not precision, past the threshold. | | Trust but verify reuse | Farm and line modes agree on all 200 map gains to $4.7\times 10^{-10}$ relative (median $7.9\times 10^{-11}$). Worst residuals over the map: $8.0\times 10^{-8}$ forward, $9.8\times 10^{-9}$ adjoint. Reuse without residuals is hope; with residuals it is a method. | | Dissipation by measurement | Leak below $10^{-11}$ (Song, Savarino meshes): $\varepsilon = 0$ is the declared choice, moving gains by at most 1.4% when raised to 0.1. Leak order one (bent meshes): dissipation required. Full procedure in the [FAQ](faq.md). | | Ranks last | 4 to 16 ranks bought 22% less time for twice the memory (18.6 to 40.6 GiB) on one point. Scale out across cases before scaling up ranks. | All figures: the linstabpy paper, computational performance section, with per-point CSVs under `examples/savarino/`. ## Reading list In the order I would read them: 1. Frantz, Loiseau and Robinet (2023), *Krylov methods for large-scale dynamical systems: application in fluid dynamics*, *Appl. Mech. Rev.* 75(3), 030802. The review behind sections 3, 5, and 6 of this page: Arnoldi factorization, Ritz convergence, non-modal bounds, and Krylov practice in fluids. 2. The [PETSc KSP manual chapter](https://petsc.org/release/manual/ksp/). Linear solves, `KSP` with `PC`, direct versus iterative, the options database. Covers sections 2 to 4 of this page with examples. 3. [`KSPPREONLY`](https://petsc.org/release/manualpages/KSP/KSPPREONLY/) and [`PCLU`](https://petsc.org/release/manualpages/PC/PCLU/). Short pages. Together they are our default solve. 4. [`MATSOLVERMUMPS`](https://petsc.org/release/manualpages/Mat/MATSOLVERMUMPS/) and the [MUMPS homepage](https://mumps-solver.org/). What the default factor backend does and its controls. 5. [`MatCreateShell`](https://petsc.org/release/manualpages/Mat/MatCreateShell/). Matrix-free operators in PETSc vocabulary. 6. [`EPSSetProblemType`](https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetProblemType/), [`EPSKRYLOVSCHUR`](https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSKRYLOVSCHUR/), [`STSINVERT`](https://slepc.upv.es/documentation/current/docs/manualpages/ST/STSINVERT/). The three pages behind every eigensolve here. 7. [petsc4py](https://petsc4py.readthedocs.io/en/stable/) and [slepc4py](https://slepc4py.readthedocs.io/en/stable/) docs. Python-side API used throughout `src/`. 8. [SciPy sparse matrices](https://docs.scipy.org/doc/scipy/reference/sparse.html). CSC and CSR formats behind the assembly side of the dual path. 9. Schmid and Henningson, *Stability and [Transition](http://www.scholarpedia.org/article/Transition_to_turbulence) in Shear Flows* (Springer, 2001). The fluids-side companion: what the eigenvalues and gains mean physically once the solvers return them. If a solver term appears elsewhere in this wiki without a link, that is a bug in that page. Link it here. --- ## References Source: https://linstabpy.vercel.app/references/ # References The papers and manuals behind this wiki, grouped by what they justify. One line each on why it is cited. To cite this software, see [how to cite](index.md#cite-linstabpy) on the home page. ## The package and its lineage * Hamada, G. Y. R. and Frantz, R. A. S. linstabpy: open-source compressible stability and resolvent analysis. *Theor. Comput. Fluid Dyn.* (in preparation, 2026). The paper for this codebase. Prefer it when citing results. * Martini, E. and Schmidt, O. T. (2024). LinStab2D: stability and resolvent analysis of compressible viscous flows in MATLAB. *Theor. Comput. Fluid Dyn.* 38, 665–685. doi:10.1007/s00162-024-00706-0. The methodology ancestor: linearized equations, SBP discretization, metrics, boundary conditions, filtering. Cite alongside the linstabpy paper. * Padovan, A., Anantharaman, V., Rowley, C. W., Vollmer, B., Colonius, T. and 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. Operator-handling design influence. ## Stability and transition theory * Schmid, P. J. and Henningson, D. S. (2001). *Stability and Transition in Shear Flows*. Springer. The textbook companion for modal spectra, branches, and transition physics. * Schmid, P. J. (2007). Nonmodal stability theory. *Annu. Rev. Fluid Mech.* 39, 129–162. Transient growth, optimal perturbations, and the finite-time viewpoint behind gain maps. * Trefethen, L. N., Trefethen, A. E., Reddy, S. C. and Driscoll, T. A. (1993). Hydrodynamic stability without eigenvalues. *Science* 261, 578–584. Why spectra alone mislead for nonnormal operators; the founding case for pseudospectra thinking. * Frantz, R. A. S., Loiseau, J.-Ch. and Robinet, J.-Ch. (2023). Krylov methods for large-scale dynamical systems: application in fluid dynamics. *Appl. Mech. Rev.* 75(3), 030802. doi:10.1115/1.4056808. Arnoldi factorization, Ritz convergence, modal vs non-modal bounds, Krylov practice. The review behind the [solver glossary](solver-glossary.md). * Hanifi, A., Schmid, P. J. and Henningson, D. S. (1996). Transient growth in compressible boundary layer flow. *Phys. Fluids* 8, 826–837. doi:10.1063/1.868864. Compressible transient growth; the wall-normal mapping parameter in the Blasius example follows this work. ## Resolvent and adjoint methods * McKeon, B. J. and Sharma, A. S. (2010). A critical-layer framework for turbulent pipe flow. *J. Fluid Mech.* 658, 336–382. Resolvent framing for amplifier flows. * Luchini, P. and Bottaro, A. (2014). Adjoint equations in stability analysis. *Annu. Rev. Fluid Mech.* 46, 493–. The adjoint machinery behind receptivity and wavemaker fields. * Chu, B.-T. (1965). On the energy transfer to small disturbances in fluid flow (Part I). *Acta Mech.* 1(3), 215–234. The compressible energy inner product used for every gain and norm here. ## Discretization * Strand, B. (1994). Summation by parts for finite difference approximations for d/dx. *J. Comput. Phys.* 110(1), 47–67. First-derivative SBP operators. * Mattsson, K. and Nordstrom, J. (2004). Summation by parts operators for finite difference approximations of second derivatives. *J. Comput. Phys.* 199(2), 503–540. Second-derivative SBP operators. * Mattsson, K., Svard, M. and Nordstrom, J. (2004). Stable and accurate artificial dissipation. *J. Sci. Comput.* 21, 57–79. doi:10.1023/B:JOMP.0000027955.75872.3f. The dissipation operator behind the epsilon knob. * Linders, V., Lundquist, T. and Nordstrom, J. (2018). On the order of accuracy of finite difference operators on diagonal norm based summation-by-parts form. *SIAM J. Numer. Anal.* 56, 1048–1063. doi:10.1137/17M1139333. Why diagonal-norm boundary closures lose order — and why only even orders 2, 4, 6, 8 ship here. * Lele, S. K. (1992). Compact finite difference schemes with spectral-like resolution. *J. Comput. Phys.* 103, 16–42. Centered stencils annihilate the sawtooth mode at any order; the reason grid-scale content needs watching. * Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. *Math. Comp.* 51, 699–706. Weight generation behind the mesh operators. ## Solver stack * Balay, S. et al. (2025). *PETSc/TAO Users Manual* (ANL-21/39, Rev. 3.24). Argonne National Laboratory. doi:10.2172/2998643. Distributed matrices, KSP/PC, the runtime under every solve. * Hernandez, V., Roman, J. E. and Vidal, V. (2005). SLEPc: a scalable and flexible toolkit for the solution of eigenvalue problems. *ACM Trans. Math. Software* 31(3), 351–362. doi:10.1145/1089014.1089019. EPS eigensolvers and spectral transformations. * Amestoy, P. R., Duff, I. S., Koster, J. and L'Excellent, J.-Y. (2001). A fully asynchronous multifrontal solver using distributed dynamic scheduling. *SIAM J. Matrix Anal. Appl.* 23(1), 15–41. MUMPS, the default sparse-direct backend. --- ## Repository README Source: https://github.com/openfluids/linstabpy#readme (External page; fetch https://github.com/openfluids/linstabpy#readme directly.)