pysbp-operators

A pure-Python (NumPy/SciPy) port of SummationByPartsOperators.jl, a Julia package by Hendrik Ranocha providing summation-by-parts (SBP) finite-difference and spectral derivative operators for provably stable discretizations of PDEs.

Status: early / partial port. This release covers periodic finite differences, Fourier pseudospectral operators, bounded-domain SBP FD, bounded and periodic upwind pairs, Lobatto--Legendre collocation, SBP-compatible dissipation, and conservative variable-coefficient diffusion. CG/DG coupling is not yet implemented. See "Roadmap" below.

No Julia runtime is required to install or use this package.

Install

pip install pysbp-operators

The command above applies after the package has been published to PyPI. Until then, install from a source checkout:

git clone https://github.com/justinh2002/pysbp-operators.git
cd pysbp-operators
pip install .

For development, including the test and quality-checking tools, use pip install -e ".[dev]".

Quickstart

import numpy as np
import pysbp_operators as sbp

D = sbp.periodic_derivative_operator(
    derivative_order=1, accuracy_order=4, xmin=0.0, xmax=2 * np.pi, N=64
)
x = sbp.grid(D)
du = D @ np.sin(x)  # ~= cos(x)

# Fourier (spectral) operator, and operator algebra:
Dx = sbp.fourier_derivative_operator(derivative_order=1, xmin=0.0, xmax=2 * np.pi, N=64)
D2 = Dx**2                      # second derivative via composition
A = D2 - 1.0                    # (d^2/dx^2 - I), still diagonal in Fourier space
f = D2 @ np.sin(3 * x) - np.sin(3 * x)
x_solved = A.solve(f)           # exact rational solve in Fourier space

sbp.integrate(np.ones_like(x), D)  # SBP quadrature: sum(M_ii * u_i)
D.to_dense(); D.to_sparse(); D.to_banded()
D.as_linear_operator()          # scipy.sparse.linalg.LinearOperator bridge

# Bounded-domain diagonal-norm SBP finite differences:
Db = sbp.derivative_operator(
    sbp.MattssonNordstrom2004(), derivative_order=1, accuracy_order=4,
    xmin=0.0, xmax=1.0, N=41,
)

# Select a member of a bounded upwind SBP bundle:
Du = sbp.upwind_operators(
    sbp.Mattsson2017, accuracy_order=6, xmin=0.0, xmax=1.0, N=41,
)
du_minus = Du.minus @ np.sin(Du.grid())

# Periodic biased Fornberg upwind pair:
Dup = sbp.upwind_operators(
    sbp.periodic_derivative_operator,
    accuracy_order=6, xmin=0.0, xmax=2 * np.pi, N=64,
)

# Negative-semidefinite artificial dissipation:
R = sbp.dissipation_operator(D, strength=0.1, order=4)
dissipative_term = R @ np.sin(x)

# Conservative variable-coefficient diffusion, ∂x(b(x) ∂x u):
L = sbp.var_coef_derivative_operator(Db, lambda x: 1.0 + x**2)
diffusive_term = L @ np.sin(L.grid())

# Lobatto--Legendre collocation (N nodes, polynomial degree N - 1):
Dl = sbp.legendre_derivative_operator(xmin=-1.0, xmax=1.0, N=12)
D2l = sbp.legendre_second_derivative_operator(xmin=-1.0, xmax=1.0, N=12)
spectral_derivative = Dl @ np.exp(Dl.grid())

Individual derivative operators are matrix-free: D @ u never forms an N x N matrix. Dense/sparse/banded conversions are provided for inspection and for composing with code that expects a matrix. An UpwindOperators bundle is intentionally ambiguous and cannot itself be multiplied by a vector; select D.minus, D.central, or D.plus.

What's implemented

• periodic_derivative_operator(...): central finite-difference operators on a periodic grid, using Fornberg's algorithm to generate minimal-width stencils of arbitrary order.
• fourier_derivative_operator(...): FFT-based spectral derivative operator.
• legendre_derivative_operator(...) and legendre_second_derivative_operator(...): Lobatto--Legendre collocation operators with exact quadrature weights and physical-interval mapping.
• derivative_operator(MattssonNordstrom2004(), ...): bounded-domain first derivatives with exact published coefficient tables at accuracy orders 2, 4, and 6, plus compatible second derivatives at accuracy orders 2, 4, and 6.
• list_sources() / lookup_source(name): extensible coefficient-source registry.
• upwind_operators(Mattsson2017, ...): bounded upwind SBP first-derivative bundles exposing .minus, .central, and .plus at accuracy orders 2, 4, and 6.
• upwind_operators(periodic_derivative_operator, ...): periodic biased Fornberg pairs at even accuracy orders, with an exact circulant central average.
• dissipation_operator(D, ...): matrix-free negative-semidefinite dissipation on every currently implemented grid family, with optional nonnegative variable strength.
• var_coef_derivative_operator(D, b): conservative matrix-free approximation of ∂x(b(x) ∂x u). Upwind bundles automatically use the dual plus/minus pair.
• grid(D), mass_matrix(D), integrate(u, D) / integrate(f, x, D).
• D.to_dense(), D.to_sparse(), D.to_banded(), D.as_linear_operator().
• Operator algebra for periodic/Fourier operators: +, -, scalar *, **, and exact rational solves like (D**2 - 1.0).solve(f), all carried out as elementwise operations on the operator's Fourier symbol (no matrix ever assembled).

Correctness

The test suite checks the algebraic identities the SBP framework depends on:

• Every norm/mass matrix is diagonal and strictly positive.
• Periodic first derivatives satisfy skew-symmetry of M @ D.
• Bounded first derivatives satisfy Q + Qᵀ = B, where Q = M @ D.
• Compatible bounded second derivatives satisfy their symmetric positive-semidefinite SBP energy decomposition.
• Mattsson (2017) upwind pairs satisfy Q+ + Q−ᵀ = B, their central member is (D− + D+) / 2, and M @ (D+ - D−) is symmetric negative-semidefinite.
• Periodic Fornberg upwind pairs satisfy the corresponding periodic dual identity with B = 0.
• Dissipation operators satisfy symmetry and negative semidefiniteness of M @ R and preserve constants.
• Variable-coefficient operators satisfy the discrete diffusion energy identity, including its endpoint flux on bounded grids.
• Lobatto--Legendre operators satisfy the first-derivative SBP identity, differentiate all polynomials representable at the nodes, and use quadrature exact through degree 2N - 3.
• Published boundary closures are checked for polynomial exactness at accuracy orders 2, 4, and 6. The periodic and Fourier operators additionally have convergence/spectral-accuracy tests.

Run the suite with:

pytest

Roadmap

Current progress:

• Phase 0 — foundations: complete.
• Phase 1 — periodic and Fourier core: substantially complete; Holoborodko differentiators remain.
• Phase 2 — bounded SBP finite differences: in progress. Mattsson--Nordstrom (2004) first and compatible second derivatives are available at accuracy orders 2, 4, and 6. Order 8 and the remaining published sources are not yet ported.
• Phase 3 — dissipation, upwind, and variable coefficients: complete for all operator families currently ported. This includes bounded Mattsson (2017) orders 2, 4, and 6, periodic biased Fornberg pairs, matrix-free SBP-compatible dissipation, and conservative variable-coefficient diffusion. Additional published coefficient-table variants remain part of Phase 2 parity work.
• Phase 4 — Legendre and CG/DG coupling: in progress. Lobatto--Legendre first and second derivatives, quadrature, variable coefficients, and dissipation are available. Uniform meshes and CG/DG coupling remain.
• Phases 5–6: AD, solver integration, performance, expanded documentation, and production release remain.

License and attribution

MIT licensed. This is an independent Python reimplementation of the algorithms and coefficient tables published in and accompanying SummationByPartsOperators.jl by Hendrik Ranocha (MIT license). Please cite the original work:

H. Ranocha, "SummationByPartsOperators.jl: A Julia library of provably stable semidiscretization techniques with mimetic properties," Journal of Open Source Software, 6(64), 3454, 2021. https://doi.org/10.21105/joss.03454

See CITATION.cff and NOTICE for details.