kernelpack-python

kernelpack-python is a Python-first port of kernelpack-matlab.

It brings the main geometry, node-generation, polynomial, RBF-FD, and solver ingredients of KernelPack into a Python codebase that is efficient, inspectable, and ready for scripting, experimentation, and extension.

This is a real port, not a placeholder wrapper around the Matlab project. The goal is to keep the same overall KernelPack-style workflow while leaning into Python strengths such as NumPy/SciPy vectorization, sparse linear algebra, and testable modular code.

The current implementation also uses targeted Numba acceleration for several hot numerical kernels, so the Python port is not limited to a purely interpretive implementation.

What it includes

Geometry objects for smooth and piecewise-smooth boundaries and surfaces: EmbeddedSurface, PiecewiseSmoothEmbeddedSurface, and RBFLevelSet
Seeded fixed-radius Poisson node generation in axis-aligned boxes, with geometry-aware clipping and boundary refinement through DomainNodeGenerator
A compact DomainDescriptor for interior, boundary, ghost nodes, normals, and nearest-neighbor search structures
Shared polynomial utilities in kernelpack.poly, including Legendre-based PolynomialBasis
RBF-FD and weighted least-squares stencil and assembly classes in kernelpack.rbffd
Numba-accelerated kernels for pairwise distances, PHS kernels, Legendre tensor-product evaluation, and stencil matrix setup
Elliptic, diffusion, and PU / PU-SL solver families in kernelpack.solvers, including PoissonSolver, VariablePoissonSolver, NonlinearVariablePoissonSolver, DiffusionSolver, PUDiffusionSolver, and the current multispecies and advection-diffusion wrappers

The main packages live in:

Supported workflows

Smooth closed curves and smooth closed 3D surfaces
Open curve segments and open 3D surface patches
Piecewise-smooth planar boundaries and piecewise 3D surfaces
Fixed-radius Poisson disk sampling in any dimension
Level-set clipping of box clouds to geometry-defined domains
Standard and overlapped RBF-FD assembly
Fixed-domain Poisson solves
Fixed-domain diffusion stepping with BDF1, BDF2, and BDF3
Variable-coefficient and nonlinear variable-coefficient Poisson solves
PU diffusion and multispecies diffusion solves
PU-SL advection and PU-SL advection-diffusion / reaction stepping

Status

The current Python port includes the main numerical layers that sit underneath typical KernelPack-style workflows:

geometry fitting and level-set construction
domain node generation
polynomial bases and multi-index utilities
local RBF-FD and weighted least-squares stencils
sparse operator assembly
fixed-domain, variable-coefficient, and nonlinear Poisson solvers
diffusion, PU diffusion, and multispecies diffusion solvers
PU-SL advection and advection-diffusion solver families
Numba-backed acceleration on the hottest geometry, polynomial, and stencil setup kernels

For overlapped RBF-FD, the current default overlap load is 0.3. That has been a better speed/accuracy tradeoff than 0.5 on the representative 3D Poisson benchmarks used during optimization.

The public API follows Python naming conventions, so Matlab-style methods such as buildClosedGeometricModelPS appear here as snake_case methods like build_closed_geometric_model_ps.

Installation

From the repository root:

python -m venv .venv
.venv\Scripts\activate
python -m pip install -e .[dev]

This installs the compiled-stack dependencies used by the package, including numpy, scipy, matplotlib, pytest, and numba.

numba is used opportunistically. The accelerated kernels are part of the normal package path, and editable installs pick them up automatically.

Run the test suite with:

python -m pytest -q

Examples and checks

Current Python checks live in:

The examples below are written as direct Python snippets so you can paste them into a script, notebook, or REPL.

Quick examples

Smooth 2D geometry

import matplotlib.pyplot as plt
import numpy as np

from kernelpack.geometry import EmbeddedSurface

# Define a smooth closed planar boundary.
t = np.linspace(0.0, 2.0 * np.pi, 50, endpoint=False)
curve = np.column_stack([np.cos(t), 0.7 * np.sin(t)])

# Build the geometric model and level set.
surface = EmbeddedSurface()
surface.set_data_sites(curve)
surface.build_closed_geometric_model_ps(2, 0.05, curve.shape[0])
surface.build_level_set_from_geometric_model()

# Extract boundary samples and normals from the fitted representation.
xb = surface.get_sample_sites()
nrmls = surface.get_nrmls()

# Plot the fitted boundary samples and normals.
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(curve[:, 0], curve[:, 1], "ko", ms=3, label="input sites")
ax.plot(xb[:, 0], xb[:, 1], ".", color="tab:teal", ms=4, label="boundary samples")
step = max(1, xb.shape[0] // 40)
ax.quiver(
    xb[::step, 0],
    xb[::step, 1],
    nrmls[::step, 0],
    nrmls[::step, 1],
    angles="xy",
    scale_units="xy",
    scale=18,
    color="tab:red",
    width=0.003,
)
ax.set_aspect("equal", adjustable="box")
ax.grid(alpha=0.2)
ax.legend(frameon=False)
plt.show()

Geometry-clipped boundary-refined nodes

import matplotlib.pyplot as plt

from kernelpack.nodes import DomainNodeGenerator

# Generate an interior-plus-boundary domain from a geometry and target spacing.
generator = DomainNodeGenerator()
domain = generator.build_domain_descriptor_from_geometry(
    surface,
    0.08,
    seed=17,
    strip_count=5,
    do_outer_refinement=True,
    outer_fraction_of_h=0.5,
    outer_refinement_zone_size_as_multiple_of_h=2.0,
)

# Pull out the packed node sets.
xi = domain.get_interior_nodes()
xb = domain.get_bdry_nodes()
xg = domain.get_ghost_nodes()

# Visualize the packed node sets.
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(xi[:, 0], xi[:, 1], ".", ms=3, color="tab:blue", label="interior")
ax.plot(xb[:, 0], xb[:, 1], ".", ms=3, color="tab:teal", label="boundary")
ax.plot(xg[:, 0], xg[:, 1], ".", ms=3, color="goldenrod", label="ghost")
ax.set_aspect("equal", adjustable="box")
ax.grid(alpha=0.2)
ax.legend(frameon=False)
plt.show()

RBF-FD operator assembly

from kernelpack.rbffd import FDDiffOp, OpProperties, RBFStencil, StencilProperties

# Ask the code to choose stencil parameters from a target accuracy.
sp = StencilProperties.from_accuracy(
    operator="lap",
    convergence_order=4,
    dimension=2,
    approximation="rbf",
    tree_mode="all",
    point_set="interior_boundary",
)

# Record stencil metadata during assembly.
op = OpProperties(record_stencils=True)

# Assemble a Laplacian on the domain descriptor.
assembler = FDDiffOp(lambda: RBFStencil())
assembler.assemble_op(domain, "lap", sp, op)
L = assembler.get_op()

For overlapped assembly, switch to FDODiffOp. The default overlap load in OpProperties is now 0.3, which is a good starting point for larger 3D problems:

from kernelpack.rbffd import FDODiffOp

op = OpProperties(record_stencils=True)  # overlap_load defaults to 0.3
assembler = FDODiffOp(lambda: RBFStencil())
assembler.assemble_op(domain, "lap", sp, op)
L = assembler.get_op()

End-to-end Poisson solve with pure Neumann data

import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import numpy as np

from kernelpack.geometry import EmbeddedSurface
from kernelpack.nodes import DomainNodeGenerator
from kernelpack.solvers import PoissonSolver

# Build a smooth closed domain.
t = np.linspace(0.0, 2.0 * np.pi, 120, endpoint=False)
curve = np.column_stack([np.cos(t), np.sin(t)])

surface = EmbeddedSurface()
surface.set_data_sites(curve)
surface.build_closed_geometric_model_ps(2, 0.06, curve.shape[0])
surface.build_level_set_from_geometric_model()

# Generate interior, boundary, and ghost nodes.
generator = DomainNodeGenerator()
domain = generator.build_domain_descriptor_from_geometry(
    surface,
    0.08,
    seed=17,
    strip_count=5,
)

# Set up an RBF-FD Poisson solve on that domain.
solver = PoissonSolver(
    lap_assembler="fd",
    bc_assembler="fd",
    lap_stencil="rbf",
    bc_stencil="rbf",
)

# Use a named target order instead of a magic number.
target_order = 4
solver.init(domain, target_order)

# Manufactured pure-Neumann problem on the unit disk.
u_exact = lambda X: (X[:, 0] ** 2 + X[:, 1] ** 2) ** 2 - (X[:, 0] ** 2 + X[:, 1] ** 2) + 1.0 / 6.0
forcing = lambda X: 4.0 - 16.0 * (X[:, 0] ** 2 + X[:, 1] ** 2)
neu_coeff = lambda Xb: np.ones(Xb.shape[0])
dir_coeff = lambda Xb: np.zeros(Xb.shape[0])
bc = lambda neu_coeffs, dir_coeffs, nr, Xb: np.sum(
    np.column_stack(
        [
            4.0 * Xb[:, 0] * (Xb[:, 0] ** 2 + Xb[:, 1] ** 2) - 2.0 * Xb[:, 0],
            4.0 * Xb[:, 1] * (Xb[:, 0] ** 2 + Xb[:, 1] ** 2) - 2.0 * Xb[:, 1],
        ]
    )
    * nr,
    axis=1,
)

# Solve and align the mean for comparison.
result = solver.solve(forcing, neu_coeff, dir_coeff, bc)
x_phys = domain.get_int_bdry_nodes()
u = result["u"]
u_true = u_exact(x_phys)
u = u - np.mean(u - u_true)
err = u - u_true

# Plot the solution and its nodal error.
tri = mtri.Triangulation(x_phys[:, 0], x_phys[:, 1])
fig, axes = plt.subplots(2, 1, figsize=(6.5, 8), constrained_layout=True)
cf0 = axes[0].tricontourf(tri, u, levels=24, cmap="viridis")
axes[0].plot(domain.get_bdry_nodes()[:, 0], domain.get_bdry_nodes()[:, 1], "k.", ms=1.5, alpha=0.5)
axes[0].set_title("Poisson solution")
fig.colorbar(cf0, ax=axes[0], shrink=0.9)

cf1 = axes[1].tricontourf(tri, err, levels=24, cmap="coolwarm")
axes[1].set_title(f"Poisson error, max |e| = {np.max(np.abs(err)):.2e}")
fig.colorbar(cf1, ax=axes[1], shrink=0.9)

for ax in axes:
    ax.set_aspect("equal", adjustable="box")
    ax.grid(alpha=0.15)

plt.show()

Diffusion stepping

import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import numpy as np

from kernelpack.solvers import DiffusionSolver

# Set up a fixed-domain diffusion stepper on the same domain.
solver = DiffusionSolver(
    lap_assembler="fd",
    bc_assembler="fd",
    lap_stencil="rbf",
    bc_stencil="rbf",
)

# Choose the diffusivity and time step.
nu = 0.25
dt = 0.02

# Use a named target order instead of a magic number.
target_order = 4
solver.init(domain, target_order, dt, nu)

# Define a manufactured transient problem with Dirichlet data.
u_exact = lambda time, X: np.exp(-time) * (X[:, 0] ** 2 + X[:, 1] ** 2)
forcing = lambda nu_value, time, X: (
    -np.exp(-time) * (X[:, 0] ** 2 + X[:, 1] ** 2)
    - 4.0 * nu_value * np.exp(-time)
)
neu_coeff = lambda Xb: np.zeros(Xb.shape[0])
dir_coeff = lambda Xb: np.ones(Xb.shape[0])
bc = lambda neu_coeffs, dir_coeffs, nr, time, Xb: u_exact(time, Xb)

# March the solution to a final time.
t_final = 0.50
nsteps = int(round(t_final / dt))
x_phys = domain.get_int_bdry_nodes()

solver.set_initial_state(u_exact(0.0, x_phys))
times = [0.0]
states = [solver.current_physical_state().copy()]

for step in range(1, nsteps + 1):
    time = step * dt
    if step == 1:
        u_next = solver.bdf1_step(time, forcing, neu_coeff, dir_coeff, bc)
    elif step == 2:
        u_next = solver.bdf2_step(time, forcing, neu_coeff, dir_coeff, bc)
    else:
        u_next = solver.bdf3_step(time, forcing, neu_coeff, dir_coeff, bc)
    times.append(time)
    states.append(u_next.copy())

# Compare against the manufactured solution at the final time.
u_final = states[-1]
u_true_final = u_exact(t_final, x_phys)
max_error = np.max(np.abs(u_final - u_true_final))
err = u_final - u_true_final

# Plot the final state, final-time error, and time history of the max nodal error.
tri = mtri.Triangulation(x_phys[:, 0], x_phys[:, 1])
fig = plt.figure(figsize=(8, 8), constrained_layout=True)
axes = fig.subplot_mosaic([["solution", "error"], ["history", "history"]])

cf0 = axes["solution"].tricontourf(tri, u_final, levels=24, cmap="viridis")
axes["solution"].set_title(f"Diffusion at t = {t_final:.2f}")
fig.colorbar(cf0, ax=axes["solution"], shrink=0.85)

cf1 = axes["error"].tricontourf(tri, err, levels=24, cmap="coolwarm")
axes["error"].set_title(f"Final-time error, max |e| = {max_error:.2e}")
fig.colorbar(cf1, ax=axes["error"], shrink=0.85)

error_history = [np.max(np.abs(state - u_exact(time, x_phys))) for time, state in zip(times, states)]
axes["history"].plot(times, error_history, color="tab:blue", lw=2)
axes["history"].scatter(times, error_history, color="tab:teal", s=18)
axes["history"].set_xlabel("t")
axes["history"].set_ylabel("max nodal error")
axes["history"].set_title("Time march to final time")
axes["history"].grid(alpha=0.2)

for key in ("solution", "error"):
    axes[key].set_aspect("equal", adjustable="box")
    axes[key].grid(alpha=0.15)

plt.show()

Package tour

`kernelpack.geometry`

This package contains the geometry-facing pieces of the port:

EmbeddedSurface for smooth closed and open fitted boundaries/surfaces
PiecewiseSmoothEmbeddedSurface for segmented piecewise boundaries
RBFLevelSet for implicit geometry representation and Newton projection
shared helpers such as distance_matrix, fibonacci_sphere, project_to_best_fit_plane, and phs_kernel

These classes are intended to support the same overall workflow as the Matlab version:

fit a geometric model to input sites
sample boundary points and normals
build an implicit level set
use that level set for node clipping and solver workflows

`kernelpack.nodes`

This package handles point generation and geometry clipping:

generate_poisson_nodes_in_box
clip_points_by_geometry
bounding_box_extents
DomainNodeGenerator

The node-generation path supports:

deterministic seeded box sampling
geometry-aware clipping
outer refinement bands near boundaries
assembly of interior, boundary, and ghost nodes into a DomainDescriptor

`kernelpack.domain`

DomainDescriptor is the compact container that links node generation and operator assembly. It stores:

interior nodes
boundary nodes
ghost nodes
boundary normals
all-node and subset node clouds
KD-tree-backed nearest-neighbor structures
domain metadata such as separation radius and level sets

`kernelpack.poly`

This package provides shared polynomial support:

Jacobi/Legendre/Chebyshev recurrence helpers
tensor-product polynomial evaluation
total-degree and hyperbolic-cross multi-index builders
PolynomialBasis, including normalization around a local center and scale

This is the polynomial backbone used by both the RBF-FD and weighted least-squares stencil paths.

`kernelpack.rbffd`

This package contains the local discretization and assembly layer:

RBFStencil
WeightedLeastSquaresStencil
StencilProperties
OpProperties
FDDiffOp
FDODiffOp

Implementation note:

Several low-level kernels in geometry, poly, and rbffd are accelerated with Numba. Solver orchestration stays in ordinary Python, while repeated numeric kernels are compiled where that pays off.

Current operator support includes:

interpolation
directional gradients
Laplacians
mixed Neumann/Dirichlet boundary-condition rows

`kernelpack.solvers`

This package contains the current solver layer:

PoissonSolver
VariablePoissonSolver
NonlinearVariablePoissonSolver
DiffusionSolver
MultiSpeciesDiffusionSolver
HeterogeneousMultiSpeciesDiffusionSolver
PUDiffusionSolver
MultiSpeciesPUDiffusionSolver
HeterogeneousMultiSpeciesPUDiffusionSolver
PUSLAdvectionSolver
MultiSpeciesPUSLAdvectionSolver
PUSLFDAdvectionDiffusionSolver
PUSLFDAdvectionDiffusionReactionSolver
PUSLPUAdvectionDiffusionSolver
PUSLPUAdvectionDiffusionReactionSolver

The solver layer stays intentionally direct. It is built on top of the same DomainDescriptor and rbffd operator assembly path rather than introducing a separate abstraction stack.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

kernelpack-python

What it includes

Supported workflows

Status

Installation

Examples and checks

Quick examples

Smooth 2D geometry

Geometry-clipped boundary-refined nodes

RBF-FD operator assembly

End-to-end Poisson solve with pure Neumann data

Diffusion stepping

Package tour

`kernelpack.geometry`

`kernelpack.nodes`

`kernelpack.domain`

`kernelpack.poly`

`kernelpack.rbffd`

`kernelpack.solvers`

FilesExpand file tree

README.md

Latest commit

History

README.md

File metadata and controls

kernelpack-python

What it includes

Supported workflows

Status

Installation

Examples and checks

Quick examples

Smooth 2D geometry

Geometry-clipped boundary-refined nodes

RBF-FD operator assembly

End-to-end Poisson solve with pure Neumann data

Diffusion stepping

Package tour

kernelpack.geometry

kernelpack.nodes

kernelpack.domain

kernelpack.poly

kernelpack.rbffd

kernelpack.solvers

`kernelpack.geometry`

`kernelpack.nodes`

`kernelpack.domain`

`kernelpack.poly`

`kernelpack.rbffd`

`kernelpack.solvers`