[WIP] Mathematical Engine Prototype for Non-Uniform Boundaries (GSoC 2026)

[utkuyilmaz1903]

Overview

Following the development of the isolated WENO engine (PR #538), this
PR provides a high-precision, zero-allocation mathematical framework
for calculating one-sided finite difference weights on non-uniform boundaries.

Handling boundaries on irregular grids is traditionally unstable and
prone to truncation errors. This engine isolates the boundary
mathematics, offering a strictly $O(1)$ and robust solution for
1st and 2nd derivatives using 4-point stencils derived via
Lagrange interpolating polynomials.

Key Features & Numerical Implementations

• Numerical Hardening (Compensated Summation): Integrates Kahan
  summation logic to calculate boundary coefficients. This ensures
  bit-perfect mass conservation ($\sum c_i \approx 0$ at the $10^{-15}$
  level), maintaining stability in long-term simulations.
• Subtraction-Free Logic: Denominators are constructed directly
  from grid intervals ($h_i$) rather than coordinate differences,
  eliminating positional reconstruction noise.
• Performance: Benchmarks show a mean execution latency of ~2.0 ns with 0 bytes allocated, ensuring zero overhead in the
  PDE solver's inner loops.
• Extreme Grid Resilience: Validated against extreme boundary layer
  stretching ratios (up to $1:10^9$) without numerical breakdown.
• MMS Verification: Verified via Method of Manufactured Solutions,
  confirming $O(\Delta x^4)$ convergence for 1st derivatives and
  $O(\Delta x^{2.6})$ for 2nd derivatives.
• AD Compatibility: Fully generic implementation (T<:Real),
  strictly supporting ForwardDiff.Dual for seamless Jacobian
  tracing in SciML architectures.
• Guardrails: Integrated assertions for grid spacing to prevent
  silent propagation of physical inconsistencies.

---

Note: Keeping this as a Draft to complement the GSoC 2026
architectural discussions regarding non-uniform grid support
in MethodOfLines.jl.

[xtalax]

Looking good, but of course we need to integrate this with the rest of the package. Which is the "centre point" of the stencil, i.e what index does it calculate? looking from the left, the missing stencil is at index 2, where index 1 would be a boundary condition.

[utkuyilmaz1903]

@xtalax The centre point here is strictly the boundary node itself (index 1 looking from the left, evaluating at $x_0$). I originally designed this to generate stencils for Neumann/Robin boundary conditions. You make a great point about index 2. The underlying Lagrange logic here can easily be shifted to evaluate the derivative at $x_1$ (index 2) instead. This would generate that missing biased stencil for the first interior point. My plan is to hook this into the rule generation flow so it can dynamically provide the correct stencil for either the boundary or that first interior point.