semiflow-pde — Python bindings for SemiFlow

Status: beta (0.15.0-beta). The API is stabilising toward 1.0; minor versions may still make breaking changes. Bug reports and feedback are welcome.

SemiFlow solves evolution equations ∂ₜu = Lu (heat, diffusion, Schrödinger, graph and manifold PDEs, …) by Chernoff approximation of operator semigroups: no matrix exponentials, no linear solves, a flat memory footprint. The core is a no_std + alloc Rust library with 3 runtime dependencies (num-traits, libm, num-complex); the same engines are available from Python, JavaScript/WebAssembly and C.

This package is the Python binding (PyO3, abi3 wheels for CPython ≥ 3.10) of the Rust crate semiflow. Every class below runs the Rust engine; the GIL is released during evolve.

The method evaluates e^{tL}f by iterating an explicit step operator S(τ): (S(t/n))ⁿ f → e^{tL} f as n → ∞. Each step is an explicit stencil over a Vec-backed grid with reused scratch buffers, so steady-state evolution does not allocate.

Mathematical foundation: Theorem 6 of I. D. Remizov (2025), Vladikavkaz Math. J. 27(4), 124–135 (doi:10.46698/a3908-1212-5385-q).

Honest performance note. The primary measured advantage is memory frugality: a flat ~3 MB working set across the 1D/2D/3D families versus 50–418 MB for heavy frameworks. Wall-clock speed is not a general strength — for matched-accuracy PDE solving SemiFlow is slower than adaptive ODE solvers and spectral methods, and parallelism pays off only on large 3D grids. Two measured niches: tail-latency-sensitive pricing (Diffusion4thChernoff ≈ 41 ns p99.9 per tick) and L1-resident low-rank carriers (TtChernoff, ReverseChernoff). Treat memory and latency as measured properties of the concrete grid types, not guarantees of the traits. The 2026-10 overhaul (docs/perf/overhaul-v0_15.md) made the grid engines 2.6–39× faster bit-for-bit, the stiff Chebyshev graph action O(√(λt)) instead of O(λt) (559× on issue #16), ETDRK4 2.1× and the general Taylor action 2.1× faster; Richardson extrapolation (richardson=) raises any engine's order by levels − 1. Head-to-head comparisons quoted in the user guide predate it.

Install

pip install semiflow-pde

The import name is semiflow. Wheels are abi3 (one wheel per platform covers CPython 3.10+). To build from a source checkout instead (needs a Rust toolchain):

pip install maturin
maturin develop --release -m crates/semiflow-py/Cargo.toml

Quickstart

Solve ∂ₜu = ∂ₓₓu from u₀(x) = e^{-x²} to t = 1 and compare with the closed-form solution u(1,x) = 5^{-1/2}·e^{-x²/5}:

import numpy as np
import semiflow

n = 1000
xs = np.linspace(-10.0, 10.0, n)
state = semiflow.Heat1D(-10.0, 10.0, n, np.exp(-xs**2))
state.evolve(t=1.0, n_steps=100)   # releases the GIL while Rust runs

exact = np.exp(-xs**2 / 5.0) / np.sqrt(5.0)
err = np.max(np.abs(state.values() - exact))
assert err < 5e-4, err

Every example on this page is executed by the project's CI against the current build.

Array I/O conventions

• All real-valued state arrays are numpy.float64 (np.float64).
• Schrödinger and SchrodingerComplex1D state arrays are numpy.complex128.
• 2D state is flat float64 in row-major x-fastest order: index j*nx + i corresponds to u(x_i, y_j).
• 3D state is flat x-fastest: index k*nx*ny + j*nx + i.
• values() always returns a copy of the internal Rust state; mutations to the returned array do not affect the object.
• Inputs are validated for NaN/Inf at construction and before evolve; non-finite inputs raise SemiflowError(kind='NanInf').
• All finite-check and grid-size errors raise SemiflowError.

Error model

All semiflow-py operations raise a single exception type:

from semiflow import SemiflowError

The .kind attribute (a string) identifies the error category:

Boundary policies

All 1D/2D/3D kernels accept a keyword argument boundary:

Far-field / Dirichlet-like boundaries for finance

Users coming from classical finite-difference option pricers often look for an inhomogeneous Dirichlet far-field ("set V = S − K e^{−rτ} at S_max") and, finding only "zero", conclude the library cannot price calls. In fact boundary="linear" is the correct and recommended closure for European call pricing.

Why "linear" works for calls: The Black-Scholes / CEV solution satisfies V(S, τ) ≈ S − K e^{−rτ} as S → ∞. This asymptote is linear in S, so V_SS → 0. Linear extrapolation from the two outermost grid nodes reproduces an affine far-field exactly (to machine precision), introducing no spurious curvature or kink. Validated on Shift1D.with_arrays with S ∈ [0, 4K], n = 1025: ATM relative error ≈ 8.5e-5 with no boundary artifact.

Puts in log-price space: When a European put is priced in log-price coordinates (z = ln S), the solution satisfies u_z → 0 at both ends of a sufficiently wide domain. Either "reflect" (zero-flux Neumann) or "linear" works; "reflect" is marginally preferable because it imposes the correct zero-derivative condition exactly.

Example — Black-Scholes call via Shift1D:

import numpy as np
import semiflow as rp

# Black-Scholes PDE: dV/dt = (sigma^2 / 2) S^2 V_SS + r S V_S - r V
# Rewrite as L V = a(S) V_SS + b(S) V_S + c(S) V with:
#   a(S) = 0.5 * sigma^2 * S^2,  b(S) = r * S,  c(S) = -r

K, T, r, sigma = 100.0, 1.0, 0.05, 0.20
S_max = 4.0 * K      # wide enough that far-field is linear
n = 1025
S_nodes = np.linspace(0.0, S_max, n)

a_arr = 0.5 * sigma**2 * S_nodes**2
b_arr = r * S_nodes
c_arr = np.full(n, -r)

# European call payoff at maturity
u0 = np.maximum(S_nodes - K, 0.0)

state = rp.Shift1D.with_arrays(
    0.0, S_max, n,
    a_arr, b_arr, c_arr,
    c_norm_bound=r,  # upper bound on |c(x)|
    u0=u0,
    boundary="linear",   # <-- correct far-field closure for calls
)
state.evolve(t=T, n_steps=200)

V = state.values()
S_atm_idx = np.argmin(np.abs(S_nodes - K))
print(f"ATM call price (semiflow): {V[S_atm_idx]:.6f}")
# rel. error vs Black-Scholes closed form ≈ 8.5e-5 at n=1025, n_steps=200

Inhomogeneous Dirichlet (future work): The Rust core provides BoundaryPolicy::Dirichlet { value } (a constant ghost-node extension), but it is not yet wired into the Python boundary= string parser. If you need to pin u = g_lo at the left edge and u = g_hi at the right edge with distinct values, use "linear" as the near-exact idiom for asymptotically-linear payoffs. Support for boundary=("dirichlet", value) is tracked in issue #20.

Usage examples

1. SchrodingerComplex1D — native complex128 wavefunction

Solve i ψ_t = (−½∂²_x + V) ψ and verify unitarity:

import numpy as np
import semiflow as rp

n = 512
xs = np.linspace(-10.0, 10.0, n)
psi0 = np.exp(-xs**2 / 2.0).astype(np.complex128)  # normalised Gaussian
dx = xs[1] - xs[0]
psi0 /= np.sqrt(np.sum(np.abs(psi0)**2) * dx)       # L2-normalise

sch = rp.SchrodingerComplex1D(-10.0, 10.0, n, psi0)
norm0 = sch.norm_squared()

sch.evolve(t=0.5, n_steps=200)

psi_t = sch.values()    # complex128 ndarray
assert abs(sch.norm_squared() / norm0 - 1.0) < 1e-12, "unitarity violated"
print(f"norm ratio = {sch.norm_squared() / norm0:.15f}")

2. Manifold2D — Riemannian manifold heat kernel

Solve ∂_t u = Δ_{S²} u on the 2-sphere via MMRS 2023 Chernoff formula:

import numpy as np
import semiflow as rp

nx, ny = 32, 64
u0 = np.zeros(nx * ny, dtype=np.float64)
u0[nx * (ny // 2) + nx // 2] = 1.0  # delta-like at chart centre

sphere = rp.Manifold2D(
    0.1, np.pi - 0.1, nx,    # theta axis
    0.0, 2 * np.pi,   ny,    # phi axis
    u0,
    manifold="sphere2",
    radius=1.0,
    curvature_correction=True,  # enables R/12 correction -> order 2
)
sphere.evolve(t=0.02, n_steps=50)

u_t = sphere.values()    # float64 ndarray, length nx*ny (row-major theta-fastest)
print(f"integral ≈ {u_t.sum() * (np.pi / nx) * (2 * np.pi / ny):.4f}")

Available manifolds: "torus" (flat T²), "sphere2" (S²(r)), "hyperbolic2" (Poincaré disk H²(s)). The radius parameter sets r or s.

3. ReactionDiffusion2D — Gray–Scott patterns, or your own kinetics

Any number of species, each with its own diffusivity; built-in kinetics by name (fisher_kpp, allen_cahn, nagumo, gray_scott, fitzhugh_nagumo, brusselator, linear) or any NumPy-vectorised callable f(t, x, u) -> du, called once per Runge–Kutta stage with all nodes (x: (dim, N), u: (K, N)), never once per node:

import numpy as np
import semiflow

n = 128
# periodic: node n is node 0, so the unit torus is sampled on [0, 1)
x, y = np.meshgrid(np.linspace(0, 1, n, endpoint=False), np.linspace(0, 1, n, endpoint=False))
bump = np.exp(-((x - 0.5)**2 + (y - 0.5)**2) / 0.005)
u0 = np.stack([1.0 - 0.5 * bump, 0.25 * bump])          # shape (K=2, ny, nx)

rd = semiflow.ReactionDiffusion2D(
    0.0, 1.0 - 1.0 / n, n, 0.0, 1.0 - 1.0 / n, n, u0,
    diffusivity=[2e-5, 1e-5], reaction="gray_scott",
    params={"feed": 0.037, "kill": 0.06}, boundary="periodic",
)
rd.evolve(2000.0, n_steps=2000)            # Strang splitting, order 2
u, v = rd.values()

# Custom kinetics + a source term, with Richardson extrapolation (order 3):
def kinetics(t, xs, u):
    return np.stack([u[0] * (1 - u[0]) - u[0] * u[1] + 0.1 * np.sin(t),
                     0.5 * u[0] * u[1] - 0.2 * u[1]])

rd2 = semiflow.ReactionDiffusion1D(-10, 10, 401, np.full((2, 401), 0.5),
                                   diffusivity=[1.0, 0.1], reaction=kinetics)
rd2.evolve(5.0, n_steps=50, richardson=2)

# Stiff chemistry (Robertson): adaptive L-stable Rosenbrock reaction steps.
def robertson(t, xs, u):
    a, b, c = u
    return np.stack([-0.04 * a + 1e4 * b * c,
                     0.04 * a - 1e4 * b * c - 3e7 * b * b, 3e7 * b * b])

u0 = np.stack([np.ones(8), np.zeros(8), np.zeros(8)])
rd3 = semiflow.ReactionDiffusion1D(0, 1, 8, u0, diffusivity=[0.0, 0.0, 0.0],
                                   reaction=robertson, reaction_solver="rosenbrock",
                                   rtol=1e-5, atol=1e-12)
rd3.evolve(40.0, n_steps=1)   # RK4 needs ~25 000 substeps per half-step here (NaN at <= 21 000)
assert abs(rd3.values()[0, 0] - 0.7158270687) < 1e-4

# Cross-diffusion: a K×K diffusivity matrix D (diagonalisable, real eigenvalues
# >= 0, verified); each diffusion step diffuses the modes of D exactly.
rd4 = semiflow.ReactionDiffusion1D(0, 1, 64, np.full((2, 64), 0.5),
                                   diffusivity=[[1e-2, 5e-3], [0.0, 2e-3]],
                                   reaction="gray_scott", params={"feed": 0.04, "kill": 0.06})
rd4.evolve(1.0, n_steps=20)

Class reference

Classes are grouped by kernel family. All stateful classes expose at least evolve(t, n_steps=100) (mutates in-place, GIL released) and values() → NDArray[np.float64] (copy). See __init__.pyi for complete signatures.

1D diffusion family

Reaction–diffusion systems (ADR-0208)

Operator splitting — multi-dimensional

Schrödinger

Both support .with_potential(v_array) and .norm_squared().

Boundary-condition kernels

Time-dependent and subordinated

Geometry and hypoelliptic operators

Graph PDE

Matrix and point-eval kernels

Anisotropic multi-D

Adjoint and adaptive wrappers

Reverse-mode AD (ADR-0156)

Constructor parameters:

.value_and_grad(tau, u0, target) -> (float, float):

import numpy as np
import semiflow as rp

n_grid = 24
xs = np.linspace(-4.0, 4.0, n_grid)

rc = rp.ReverseHeat1D(theta=0.4, xmin=-4.0, xmax=4.0, n_grid=n_grid, n_steps=8)
u0     = np.exp(-xs**2)
target = np.zeros(n_grid)

value, grad = rc.value_and_grad(tau=0.05, u0=u0, target=target)
print(f"loss={value:.6e}  ∂J/∂θ={grad:.6e}")

Raises SemiflowError with .kind in {'OutOfDomain', 'GridMismatch', 'NanInf'}.

NARROW scope (§51.5): constant-a DiffusionChernoff only; θ is the uniform diffusivity. Variable-coefficient and nonlinear kernels are out of scope for ReverseHeat1D.

Tensor-train carriers (ADR-0171 / ADR-0178)

Curse-escaped O(d·n·r²) storage for separable diagonal-A diffusion on ℝᵈ. State and evolvers share the TtState carrier; TtEvolver / TtCoupledEvolver cover constant-coefficient cases, VarCoefTtEvolver covers variable-coefficient additive-separable operators (issue #2, ADR-0178).

Gridless / particle carriers (ADR-0171)

Sparse signed-measure ensemble on ℝ (D=1). Curse-escaped: the 3ᴰ dense tree is never materialised; only sparse marginals and scalar observables cross the Python boundary.

v3 Evolver surface

Greeks and hyper-dual AD

Adjoint Fokker-Planck (particle pushforward)

Killing and absorbing boundary conditions

Obstacle / variational-inequality family

Resolvent time-jump family (TWS contour)

Matrix diffusion (coupled 2-component)

Wentzell dynamic boundary conditions

Additional 1D kernels

Sparse-grid and high-dimensional kernels

Graph adjoint and Krylov families

Symmetric operator and FEM utilities

Example — stiff operator: the rational and implicit paths. The polynomial Krylov paths need O(√(λ_max·t)) mat-vecs; path="rational" (shift-and-invert Lanczos, ADR-0206) reaches tol in a fixed number of O(N) tridiagonal solves (27 at tol = 1e-10) whatever the stiffness. path="implicit" (Jacobi-preconditioned backward-Euler PCG, accuracy O(t/n_steps), about √κ CG iterations per sub-step with κ ≈ 1 + (t/n_steps)·λ_max) handles stiffness such as λ_max ≈ 4·10⁷ for smooth data. Rough data (e.g. white noise) at this stiffness needs more CG iterations than the fixed cap of 400 and raises SemiflowError with kind == "ConvergenceFailed".

import numpy as np
import semiflow

# Stiff 1-D Neumann Laplacian (tridiag(-1, 2, -1), boundary rows 1, -1) * 1e7, N = 400.
n, scale = 400, 1e7
indptr, indices, data = [0], [], []
for i in range(n):
    for j in (i - 1, i, i + 1):
        if 0 <= j < n:
            indices.append(j)
            data.append(scale * ((1.0 if i in (0, n - 1) else 2.0) if j == i else -1.0))
    indptr.append(len(indices))

op = semiflow.SymmetricOperator.from_csr(
    np.array(indptr, dtype=np.int64), np.array(indices, dtype=np.int32),
    np.array(data), n, 1e-6)
v0 = np.linspace(1.0 / n, 1.0, n).reshape(n, 1)
u = op.evolve_batched(1.0, v0, path="implicit", n_steps=100)

# Everything but the constant (null-space) mode has decayed: u -> mean(v0).
assert np.max(np.abs(u - v0.mean())) < 1e-6

# The rational path: 27 tridiagonal solves, error below 1e-10.
u = op.evolve_batched(1.0, v0, path="rational")
assert np.max(np.abs(u - v0.mean())) < 1e-10

mass_lumped_evolve accepts the same path= keyword ("rational" uses M directly); MassKOperator.evolve the polynomial and implicit ones.

Free functions

version(), heisenberg_heat_kernel(), sample_gridfn2d(), edge_weight_grad(), edge_weight_grad_batched(), graph_expmv_frechet(), symmetric_op_expmv_frechet(), assemble_conservative_csr_1d(), mass_lumped_evolve(), phi_action(), phi_action_batched(), phi_combination(), shift1d_coeff_grad(), richardson_weights() — see __init__.pyi for their signatures.

Performance

GIL release follows the three-phase py.detach pattern (ADR-0031): acquire → snapshot inputs → detach → Rust compute → reacquire. Send + Sync is verified at compile time with static_assertions.

Indicative timings on i7-12700K (1000 nodes, 100 steps, Heat1D):

.with_a_array uses a pure-Rust Arc<Vec<f64>> Catmull-Rom interpolant and never re-acquires the GIL during evolve. .with_a_function calls the Python callables about 10 times per node (a twice, a_prime 7×, a_double_prime once) per τ-keyed step plan (ADR-0215). There is one plan per evolve call of n ≥ 2 steps, per chunk of evolve_chunked (of ≥ 2 steps) and per Richardson level, not one per step. Single-step calls, and grids whose plan exceeds 2²² doubles, pay the callbacks on every step; prefer .with_a_array for those.

Type stubs

__init__.pyi and the py.typed marker ship with every wheel. Static type checkers (mypy, pyright, pylance) pick them up automatically. The pyrightconfig.json at the repo root adds crates/semiflow-py/python to extraPaths so local development also resolves the stubs correctly (0 reportAttributeAccessIssue errors).

Bindings

The bindings mirror the user-facing engines, not the internal composition types (AxisLift, StrangSplit, … are Rust-side building blocks). Variable coefficients cross the language boundary as pre-sampled arrays (evaluated in Rust, no call back into the host language during evolve); some classes also accept a host-language callback (e.g. Heat1D.with_a_function in Python, Heat1D.withAFunction in JS, smf_state_new_with_closure in C). The Bindings guide has side-by-side examples.

Documentation

Accuracy and design

Every numerical claim is gated in CI against a closed-form or high-order reference oracle (convergence-order and sup-norm tests). Design principles: no_std + alloc core; 3 runtime dependencies; SIMD hot paths isolated to src/simd/ (AVX2 on x86_64, NEON on aarch64, scalar fallback elsewhere); bit-reproducible parallelism across thread counts; no unsafe in the math kernels. Builds with and without simd agree to rounding but not always bit for bit, so golden files are tied to one configuration.

How to cite

If you use SemiFlow in academic work, cite both the software and the theorem.

@software{volkov2026semiflow,
  author  = {Volkov, Ilia},
  title   = {{SemiFlow}: {Chernoff} Approximation of Operator Semigroups},
  year    = {2026},
  version = {0.15.0-beta},
  doi     = {10.5281/zenodo.20837851},
  url     = {https://doi.org/10.5281/zenodo.20837851}
}

@article{Remizov2025,
  author  = {I. D. Remizov},
  title   = {Chernoff Approximations of the Solution of Linear ODE with Variable Coefficients},
  journal = {Vladikavkaz Math. J.},
  volume  = {27},
  number  = {4},
  pages   = {124--135},
  year    = {2025},
  doi     = {10.46698/a3908-1212-5385-q}
}

The concept DOI 10.5281/zenodo.20837851 always resolves to the latest version; CITATION.cff has the full entry (GitHub's "Cite this repository" reads it).

License

Dual-licensed under MIT or Apache-2.0, at your option.

Author, contributing, security

SemiFlow is created and maintained by Ilia Volkov. See CONTRIBUTING.md and CODE_OF_CONDUCT.md. Report vulnerabilities privately as described in SECURITY.md — not in public issues.