PDE Layer (pypielm.pde)

PDE building blocks: operators, collocation, and constraints.

Public surface:

from pypielm.pde import (
    # Operators
    gradient, laplacian, divergence, advection_term, AnalyticLaplacian,
    # Domain + samplers
    BoxDomain, UniformSampler, LHSSampler, AdaptiveSampler, GridSampler,
    # Conditions
    DirichletBC, NeumannBC, InitialCondition, PeriodicBC,
)
pypielm.pde.gradient(u, x)[source]

Compute the gradient ∇u w.r.t. x via autograd.

\[\nabla u = \left(\frac{\partial u}{\partial x_1}, \ldots, \frac{\partial u}{\partial x_d}\right)\]
Parameters:

  • u (Tensor) – Scalar field tensor, shape (N,) or (N, 1). Must have requires_grad=True set on the computation graph leading to it.

  • x (Tensor) – Input coordinates, shape (N, d), with requires_grad=True.

Return type:

Tensor

Returns:

Gradient tensor of shape (N, d).

pypielm.pde.laplacian(u, x)[source]

Compute the Laplacian Δu = Σᵢ ∂²u/∂xᵢ² via autograd.

\[\Delta u = \sum_{i=1}^{d} \frac{\partial^2 u}{\partial x_i^2}\]
Parameters:

  • u (Tensor) – Scalar field, shape (N,) or (N, 1).

  • x (Tensor) – Input coordinates, shape (N, d), with requires_grad=True.

Return type:

Tensor

Returns:

Laplacian values, shape (N, 1).

pypielm.pde.divergence(flux, x)[source]

Compute the divergence ∇ · F of a vector field F.

\[\nabla \cdot \mathbf{F} = \sum_{i=1}^{d} \frac{\partial F_i}{\partial x_i}\]
Parameters:

  • flux (Tensor) – Vector field tensor, shape (N, d).

  • x (Tensor) – Input coordinates, shape (N, d), with requires_grad=True.

Return type:

Tensor

Returns:

Divergence values, shape (N, 1).

pypielm.pde.advection_term(u, v, x)[source]

Compute the advection term v · ∇u.

\[\mathbf{v} \cdot \nabla u = \sum_{i=1}^{d} v_i \frac{\partial u}{\partial x_i}\]
Parameters:

  • u (Tensor) – Scalar field, shape (N,) or (N, 1).

  • v (Tensor) – Advection velocity, shape (N, d).

  • x (Tensor) – Input coordinates, shape (N, d), with requires_grad=True.

Return type:

Tensor

Returns:

Advection term, shape (N, 1).

class pypielm.pde.AnalyticLaplacian(feature_map=None, input_dim=None)[source]

Fast Laplacian operator using precomputed analytic second derivatives.

Avoids the O(N · H · d) autograd overhead by using the analytic d2() (or laplacian()) method of the feature map directly.

When feature_map is provided and exposes a laplacian method, that is used directly. Otherwise the d2() method is summed over all spatial dimensions.

\[\Delta [H \boldsymbol{\beta}](x) = \left(\sum_{i=1}^{d} \mathbf{H}^{(2,i)}(x)\right) \boldsymbol{\beta}\]

where \(\mathbf{H}^{(2,i)}\) is the second partial derivative of H w.r.t. coordinate \(x_i\).

Parameters:

  • feature_map (RandomFeatureMap | None) – A RandomFeatureMap or compatible feature map instance (must expose d2 or laplacian).

  • input_dim (int | None) – Number of spatial dimensions to sum over. Inferred from feature_map.input_dim if feature_map is provided.

class pypielm.pde.BoxDomain(lb, ub)[source]

Axis-aligned bounding box in R^d.

Parameters:

Example:

domain = BoxDomain(lb=[0.0, 0.0], ub=[1.0, 1.0])
property dim: int

Spatial dimension d.

class pypielm.pde.UniformSampler(domain, n_points=1000, seed=42)[source]

Sample collocation points uniformly at random within a BoxDomain.

Parameters:

  • domain (BoxDomain) – The spatial domain.

  • n_points (int) – Number of collocation points to sample.

  • seed (int) – Random seed.

sample()[source]

Draw n_points uniform samples from domain.

Return type:

Tensor

Returns:

Tensor of shape (n_points, d).

class pypielm.pde.LHSSampler(domain, n_points=1000, seed=42)[source]

Latin Hypercube Sampling (LHS) within a BoxDomain.

LHS ensures better space-filling than pure uniform sampling. Uses scipy.stats.qmc.LatinHypercube when available; falls back to a stratified uniform sampler otherwise.

Parameters:

  • domain (BoxDomain) – The spatial domain.

  • n_points (int) – Number of collocation points.

  • seed (int) – Random seed.

sample()[source]

Draw n_points LHS samples.

Return type:

Tensor

Returns:

Tensor of shape (n_points, d).

class pypielm.pde.AdaptiveSampler(domain, residual_fn, n_points=1000, refine_ratio=0.5, n_candidates=10000, seed=42)[source]

Residual-guided adaptive collocation sampler.

Samples a large candidate set, evaluates the provided residual_fn, and returns points concentrated in high-residual regions (top refine_ratio fraction by residual magnitude) padded with uniform samples.

Parameters:

  • domain (BoxDomain) – The spatial domain.

  • residual_fn (Callable[[Tensor], Tensor]) – Callable f(X: Tensor) → Tensor returning scalar residuals of shape (N,) or (N, 1) for input of shape (N, d).

  • n_points (int) – Number of collocation points to return.

  • refine_ratio (float) – Fraction of returned points that are residual-guided (the rest are uniform samples).

  • n_candidates (int) – Number of candidates to evaluate before selecting.

  • seed (int) – Random seed.

sample()[source]

Sample collocation points with residual-guided refinement.

Return type:

Tensor

Returns:

Tensor of shape (n_points, d).

class pypielm.pde.GridSampler(domain, nx=64, ny=64)[source]

Structured Cartesian grid sampler (1D, 2D, or higher).

Useful for finite-difference baselines and structured visualisation.

Parameters:

  • domain (BoxDomain) – The spatial domain.

  • nx (int) – Number of grid points along axis 0.

  • ny (int) – Number of grid points along axis 1 (ignored for 1D domains).

sample()[source]

Return all grid points as a tensor.

Return type:

Tensor

Returns:

Tensor of shape (N_total, d) where N_total = nx (1D) or nx * ny (2D).

class pypielm.pde.DirichletBC(boundary_fn, points, weight=1.0)[source]

Hard Dirichlet boundary condition: u(x) = g(x) on ∂Ω.

Parameters:

  • boundary_fn (Callable[[Tensor], Tensor]) – Callable g(x: Tensor) → Tensor returning the prescribed values at boundary points, shape (N_bc,) or (N_bc, 1).

  • points (Tensor) – Boundary collocation points, shape (N_bc, d).

  • weight (float) – Observation precision for this BC block.

assemble(feature_map)[source]

Evaluate BC and return the linear system block.

Parameters:

feature_map (RandomFeatureMap) – The model’s hidden-layer feature map.

Return type:

WeightedLinearSystem

Returns:

WeightedLinearSystem with H = feature_map(points), y = boundary_fn(points), weight = self.weight.

class pypielm.pde.NeumannBC(flux_fn, normal, points, weight=1.0)[source]

Neumann boundary condition: ∂u/∂n = h(x) on ∂Ω.

The feature-matrix contribution uses the outward unit normal n to form the directional derivative: ∂H/∂n = Σᵢ nᵢ · (∂H/∂xᵢ).

Parameters:

  • flux_fn (Callable[[Tensor], Tensor]) – Callable returning the prescribed normal flux, shape (N_bc,) or (N_bc, 1).

  • normal (Tensor) – Outward unit normal vectors, shape (N_bc, d).

  • points (Tensor) – Boundary collocation points, shape (N_bc, d).

  • weight (float) – Observation precision.

assemble(feature_map)[source]

Evaluate flux BC and return the linear system block.

Builds the directional-derivative feature matrix:

H_n[i, j] = Σ_k  n[i, k] * (∂H/∂x_k)[i, j]
Return type:

WeightedLinearSystem

class pypielm.pde.InitialCondition(ic_fn, points, weight=1.0)[source]

Initial condition: u(x, t=0) = u₀(x).

Parameters:

  • ic_fn (Callable[[Tensor], Tensor]) – Callable returning prescribed initial values, shape (N_ic,) or (N_ic, 1).

  • points (Tensor) – Initial condition points (with t=0 already embedded), shape (N_ic, d).

  • weight (float) – Observation precision.

assemble(feature_map)[source]

Evaluate IC and return the linear system block.

Return type:

WeightedLinearSystem

class pypielm.pde.PeriodicBC(axis, points_left, points_right, weight=1.0)[source]

Periodic boundary condition along a specified axis.

Pairs boundary points x_left and x_right and enforces u(x_left) = u(x_right) by adding the penalty rows H(x_left) - H(x_right) with target y = 0 to the linear system.

Parameters:

  • axis (int) – Axis along which periodicity is imposed (informational only; the caller is responsible for pairing points correctly).

  • points_left (Tensor) – Left boundary points, shape (N_bc, d).

  • points_right (Tensor) – Right boundary points, shape (N_bc, d).

  • weight (float) – Observation precision for the pairing rows.

assemble(feature_map)[source]

Assemble the pairing penalty block.

Returns rows H(x_left) - H(x_right) with target zero, so the solver enforces u(x_left) ≈ u(x_right).

Return type:

WeightedLinearSystem

Differential Operators

PDE differential operators (autograd + analytic fast paths).

Public API:

from pypielm.pde.operators import (
    gradient, laplacian, divergence, advection_term, AnalyticLaplacian
)
pypielm.pde.operators.gradient(u, x)[source]

Compute the gradient ∇u w.r.t. x via autograd.

\[\nabla u = \left(\frac{\partial u}{\partial x_1}, \ldots, \frac{\partial u}{\partial x_d}\right)\]
Parameters:

  • u (Tensor) – Scalar field tensor, shape (N,) or (N, 1). Must have requires_grad=True set on the computation graph leading to it.

  • x (Tensor) – Input coordinates, shape (N, d), with requires_grad=True.

Return type:

Tensor

Returns:

Gradient tensor of shape (N, d).

pypielm.pde.operators.laplacian(u, x)[source]

Compute the Laplacian Δu = Σᵢ ∂²u/∂xᵢ² via autograd.

\[\Delta u = \sum_{i=1}^{d} \frac{\partial^2 u}{\partial x_i^2}\]
Parameters:

  • u (Tensor) – Scalar field, shape (N,) or (N, 1).

  • x (Tensor) – Input coordinates, shape (N, d), with requires_grad=True.

Return type:

Tensor

Returns:

Laplacian values, shape (N, 1).

pypielm.pde.operators.divergence(flux, x)[source]

Compute the divergence ∇ · F of a vector field F.

\[\nabla \cdot \mathbf{F} = \sum_{i=1}^{d} \frac{\partial F_i}{\partial x_i}\]
Parameters:

  • flux (Tensor) – Vector field tensor, shape (N, d).

  • x (Tensor) – Input coordinates, shape (N, d), with requires_grad=True.

Return type:

Tensor

Returns:

Divergence values, shape (N, 1).

pypielm.pde.operators.advection_term(u, v, x)[source]

Compute the advection term v · ∇u.

\[\mathbf{v} \cdot \nabla u = \sum_{i=1}^{d} v_i \frac{\partial u}{\partial x_i}\]
Parameters:

  • u (Tensor) – Scalar field, shape (N,) or (N, 1).

  • v (Tensor) – Advection velocity, shape (N, d).

  • x (Tensor) – Input coordinates, shape (N, d), with requires_grad=True.

Return type:

Tensor

Returns:

Advection term, shape (N, 1).

class pypielm.pde.operators.AnalyticLaplacian(feature_map=None, input_dim=None)[source]

Bases: object

Fast Laplacian operator using precomputed analytic second derivatives.

Avoids the O(N · H · d) autograd overhead by using the analytic d2() (or laplacian()) method of the feature map directly.

When feature_map is provided and exposes a laplacian method, that is used directly. Otherwise the d2() method is summed over all spatial dimensions.

\[\Delta [H \boldsymbol{\beta}](x) = \left(\sum_{i=1}^{d} \mathbf{H}^{(2,i)}(x)\right) \boldsymbol{\beta}\]

where \(\mathbf{H}^{(2,i)}\) is the second partial derivative of H w.r.t. coordinate \(x_i\).

Parameters:

  • feature_map (RandomFeatureMap | None) – A RandomFeatureMap or compatible feature map instance (must expose d2 or laplacian).

  • input_dim (int | None) – Number of spatial dimensions to sum over. Inferred from feature_map.input_dim if feature_map is provided.

Collocation Samplers

Collocation point samplers and domain descriptors.

Public API:

from pypielm.pde.collocation import (
    BoxDomain, UniformSampler, LHSSampler, AdaptiveSampler, GridSampler
)
class pypielm.pde.collocation.BoxDomain(lb, ub)[source]

Bases: object

Axis-aligned bounding box in R^d.

Parameters:

Example:

domain = BoxDomain(lb=[0.0, 0.0], ub=[1.0, 1.0])
property dim: int

Spatial dimension d.

class pypielm.pde.collocation.UnionDomain(domains)[source]

Bases: object

Union of multiple BoxDomain objects.

Parameters:

domains (list[BoxDomain]) – List of BoxDomain objects forming the union.

property dim: int
class pypielm.pde.collocation.UniformSampler(domain, n_points=1000, seed=42)[source]

Bases: object

Sample collocation points uniformly at random within a BoxDomain.

Parameters:

  • domain (BoxDomain) – The spatial domain.

  • n_points (int) – Number of collocation points to sample.

  • seed (int) – Random seed.

sample()[source]

Draw n_points uniform samples from domain.

Return type:

Tensor

Returns:

Tensor of shape (n_points, d).

class pypielm.pde.collocation.LHSSampler(domain, n_points=1000, seed=42)[source]

Bases: object

Latin Hypercube Sampling (LHS) within a BoxDomain.

LHS ensures better space-filling than pure uniform sampling. Uses scipy.stats.qmc.LatinHypercube when available; falls back to a stratified uniform sampler otherwise.

Parameters:

  • domain (BoxDomain) – The spatial domain.

  • n_points (int) – Number of collocation points.

  • seed (int) – Random seed.

sample()[source]

Draw n_points LHS samples.

Return type:

Tensor

Returns:

Tensor of shape (n_points, d).

class pypielm.pde.collocation.AdaptiveSampler(domain, residual_fn, n_points=1000, refine_ratio=0.5, n_candidates=10000, seed=42)[source]

Bases: object

Residual-guided adaptive collocation sampler.

Samples a large candidate set, evaluates the provided residual_fn, and returns points concentrated in high-residual regions (top refine_ratio fraction by residual magnitude) padded with uniform samples.

Parameters:

  • domain (BoxDomain) – The spatial domain.

  • residual_fn (Callable[[Tensor], Tensor]) – Callable f(X: Tensor) → Tensor returning scalar residuals of shape (N,) or (N, 1) for input of shape (N, d).

  • n_points (int) – Number of collocation points to return.

  • refine_ratio (float) – Fraction of returned points that are residual-guided (the rest are uniform samples).

  • n_candidates (int) – Number of candidates to evaluate before selecting.

  • seed (int) – Random seed.

sample()[source]

Sample collocation points with residual-guided refinement.

Return type:

Tensor

Returns:

Tensor of shape (n_points, d).

class pypielm.pde.collocation.GridSampler(domain, nx=64, ny=64)[source]

Bases: object

Structured Cartesian grid sampler (1D, 2D, or higher).

Useful for finite-difference baselines and structured visualisation.

Parameters:

  • domain (BoxDomain) – The spatial domain.

  • nx (int) – Number of grid points along axis 0.

  • ny (int) – Number of grid points along axis 1 (ignored for 1D domains).

sample()[source]

Return all grid points as a tensor.

Return type:

Tensor

Returns:

Tensor of shape (N_total, d) where N_total = nx (1D) or nx * ny (2D).

Boundary & Initial Conditions

Boundary and initial condition helpers.

Each condition class evaluates the constraint at given boundary/IC points and returns a WeightedLinearSystem tuple (H_bc, y_bc, weight) that can be stacked into the global linear system assembled during model training.

Public API:

from pypielm.pde.constraints import (
    DirichletBC, NeumannBC, InitialCondition, PeriodicBC
)
class pypielm.pde.constraints.DirichletBC(boundary_fn, points, weight=1.0)[source]

Bases: object

Hard Dirichlet boundary condition: u(x) = g(x) on ∂Ω.

Parameters:

  • boundary_fn (Callable[[Tensor], Tensor]) – Callable g(x: Tensor) → Tensor returning the prescribed values at boundary points, shape (N_bc,) or (N_bc, 1).

  • points (Tensor) – Boundary collocation points, shape (N_bc, d).

  • weight (float) – Observation precision for this BC block.

assemble(feature_map)[source]

Evaluate BC and return the linear system block.

Parameters:

feature_map (RandomFeatureMap) – The model’s hidden-layer feature map.

Return type:

WeightedLinearSystem

Returns:

WeightedLinearSystem with H = feature_map(points), y = boundary_fn(points), weight = self.weight.

class pypielm.pde.constraints.NeumannBC(flux_fn, normal, points, weight=1.0)[source]

Bases: object

Neumann boundary condition: ∂u/∂n = h(x) on ∂Ω.

The feature-matrix contribution uses the outward unit normal n to form the directional derivative: ∂H/∂n = Σᵢ nᵢ · (∂H/∂xᵢ).

Parameters:

  • flux_fn (Callable[[Tensor], Tensor]) – Callable returning the prescribed normal flux, shape (N_bc,) or (N_bc, 1).

  • normal (Tensor) – Outward unit normal vectors, shape (N_bc, d).

  • points (Tensor) – Boundary collocation points, shape (N_bc, d).

  • weight (float) – Observation precision.

assemble(feature_map)[source]

Evaluate flux BC and return the linear system block.

Builds the directional-derivative feature matrix:

H_n[i, j] = Σ_k  n[i, k] * (∂H/∂x_k)[i, j]
Return type:

WeightedLinearSystem

class pypielm.pde.constraints.InitialCondition(ic_fn, points, weight=1.0)[source]

Bases: object

Initial condition: u(x, t=0) = u₀(x).

Parameters:

  • ic_fn (Callable[[Tensor], Tensor]) – Callable returning prescribed initial values, shape (N_ic,) or (N_ic, 1).

  • points (Tensor) – Initial condition points (with t=0 already embedded), shape (N_ic, d).

  • weight (float) – Observation precision.

assemble(feature_map)[source]

Evaluate IC and return the linear system block.

Return type:

WeightedLinearSystem

class pypielm.pde.constraints.PeriodicBC(axis, points_left, points_right, weight=1.0)[source]

Bases: object

Periodic boundary condition along a specified axis.

Pairs boundary points x_left and x_right and enforces u(x_left) = u(x_right) by adding the penalty rows H(x_left) - H(x_right) with target y = 0 to the linear system.

Parameters:

  • axis (int) – Axis along which periodicity is imposed (informational only; the caller is responsible for pairing points correctly).

  • points_left (Tensor) – Left boundary points, shape (N_bc, d).

  • points_right (Tensor) – Right boundary points, shape (N_bc, d).

  • weight (float) – Observation precision for the pairing rows.

assemble(feature_map)[source]

Assemble the pairing penalty block.

Returns rows H(x_left) - H(x_right) with target zero, so the solver enforces u(x_left) ≈ u(x_right).

Return type:

WeightedLinearSystem