Kernels and Structures

The following sections offer a basic overview of the mathematical background for some of the core concepts in smudgy. First, we introduce the theory of moment tensors, that are estimated from a discrete sample of points. Then, we introduce the concept of kernels, how they are defined in smudgy and how their shape relates to the local estimates of the moment tensors.

Kernel structures and moment tensors

Estimating continuous distributions from discrete particle samples is a central problem in particle-based methods such as SPH. A common approach is to locally approximate the underlying distribution around each particle by defining a smoothing kernel and an associated distance measure \(d(x)\).

Suppose we have a point cloud \(C\) containing particles \(p \in C\) with coordinates \(x_p\) and weights \(m_p\). \(C\) could for instance be the result of a nieghbor search operation. The total mass of \(C\) is denoted as \(M = \sum_p m_p\). For any given point \(\mu\), we can compute the local geometry by using \(C\) and estimating the true underlying distribution of particles. We can do this with a kernel and the distance measure \(d(x)\):

• separable: In the simplest case, a 1D kernel function is defined for each dimension and the total kernel is simply the product of all its 1D shapes, typically resulting in rectangular shaped kernels.

• isotropic: The next higher order assumes isotropic support regions, where the kernel depends only on the distance to the particle center. Thus, for the case of separable and isotropic kernel shapes, the kernel distance measure is defined via the Eucledian distance as

  \[ d(x) = ||x - \mu||, \]

  where \(\mu\) is the local particle center. This definition leads to circular (2D) or spherical (3D) kernel shapes.

• anisotropic: The distance measure in the second-order method is the Mahalanobis distance, defined via the covariance tensor \(\Sigma\):

  \[ d^2(x) = (x-\mu)^T \Sigma^{-1} (x-\mu) \qquad \text{where} \qquad \Sigma_{ij} = \frac{1}{M} \sum_p m_p (x_{p,i}-\mu_i) (x_{p,j}-\mu_j), \]

  The resulting kernel support and shape becomes elliptical (2D) or ellipsoidal (3D).

  Caution

  The current implementation estimates covariance tensors from a standard radial neighbor search and does not iteratively refine the neighbor set using the anisotropic metric.

• higher orders: Higher order methods estimate local geometric structure from neighboring particles, allowing anisotropic and potentially higher-order kernel shapes. In general, the implicit kernel distance may be written as

  \[ d^2(x) = \sum_{n=2}^{N} \sum_{i_1,\dots,i_n=1}^{d} M^{(n)}_{i_1\dots i_n} \prod_{k=1}^{n} (x_{i_k}-\mu_{i_k}) \]

  where the \(n\)-th order moment tensor \(M^{(n)}\) is estimated from neighboring particles via:

  \[ M^{(n)}_{i_1 \dots i_n} = \frac{1}{M} \sum_p m_p \prod_{k=1}^{n} (x_{p,i_k} - \mu_{i_k}). \]

  For instance, the third-order tensor (skewness) introduces asymmetric distortions,

  \[ M^{(3)} = S_{ijk} = \frac{1}{M} \sum_p m_p (x_{p,i}-\mu_i) (x_{p,j}-\mu_j) (x_{p,k}-\mu_k), \]

  while the fourth-order tensor (kurtosis) captures higher-order shape structure such as tail behavior and nonlinear deformations,

  \[ M^{(4)}= K_{ijkl} = \frac{1}{M} \sum_p m_p (x_{p,i}-\mu_i) (x_{p,j}-\mu_j) (x_{p,k}-\mu_k) (x_{p,l}-\mu_l). \]

  These higher-order terms then affect the distance measure via

  \[ d^2(x) = (x-\mu)^T \Sigma^{-1} (x-\mu) + S_{ijk}\Delta x_i \Delta x_j \Delta x_k + K_{ijkl}\Delta x_i \Delta x_j \Delta x_k \Delta x_l +\dots \]

  with \(\Delta x_i = x_i - \mu_i\).

Note

Currently, the package supports separable, isotropic and anisotropic covariance-based kernels, while higher-order moment methods are planned for future releases.

Kernels in smudgy

In smudgy, kernels \(W\) are defined via the smoothing lengths or tensors, a normalization constant \(\sigma\) (typicall dimension dependent) and a shape function \(K(q)\). Their evaluation, according to the distance measures introduced above, depend on the chosen strucutre:

• separable kernels are defined through their 1D counterparts. This option enables analytically exact integrals when depositing particle fields onto grids. With the spatial dimension \(d\) they are defined as

\[ W(\mathbf r; h) = \prod_{i=1}^{d} W_{\rm 1D}(x_i; h), \]

• isotropic kernels with smoothing length \(h\), \(\sigma_d\) a dimension‑dependent normalization and \(K(q)\) the shape function with normalized coordinates \(q\) are defined as

\[ W(\mathbf r; h) = \frac{\sigma_d}{h^{d}}\,K\!\left(q\right), \qquad q=\frac{\|\mathbf r\|}{h}, \]

• anisotropic kernels are defined via the smoothing matrix \(H = \sqrt{\Sigma}\), so that the above expression generalizes to

\[ W(\mathbf r; H) = \frac{\sigma_d}{\det(H)}\,K\big( \xi \big), \qquad \xi = \|H^{-1}\mathbf r\| \]

The table below summarizes the available kernel types in smudgy:

Name	string	\(K(q)\) shape	support	\(\sigma_\mathrm{1D}\)	\(\sigma_\mathrm{2D}\)	\(\sigma_\mathrm{3D}\)
Tophat	'tophat'	\(\begin{cases}1 & q \le 0.5 \\ 0 & q > 0.5\end{cases}\)	\(0.5\)	\(1\)	\(4/\pi\)	\(6/\pi\)
TSC	'tsc'	\(\begin{cases}\frac{3}{4}-q^2 & q < 0.5 \\ \frac{1}{2}(1.5-q)^2 & 0.5 \le q < 1.5 \\ 0 & q \ge 1.5\end{cases}\)	\(1.5\)	\(1\)	\(3/\pi\)	\(3/\pi\)
Gaussian	'gaussian'	\(\exp(-q^2)\)	\(3.0\)	\(1/\sqrt{\pi}\)	\(1/\pi\)	\(1/\pi^{3/2}\)
Lucy	'lucy'	\(\begin{cases}(1+3q)(1-q)^3 & q \le 1 \\ 0 & q > 1\end{cases}\)	\(1.0\)	\(5/4\)	\(5/\pi\)	\(105/(16\pi)\)
Cubic Spline	'cubic_spline'	\(\begin{cases}1 - 6q^2 + 6q^3 & 0 \le q < 0.5 \\ 2(1-q)^3 & 0.5 \le q < 1 \\ 0 & q \ge 1\end{cases}\)	\(1.0\)	\(4/3\)	\(40/(7\pi)\)	\(8/\pi\)
Quintic Spline	'quintic_spline'	\(\begin{cases}(3-q)^5 - 6(2-q)^5 + 15(1-q)^5 & 0 \le q < 1 \\ (3-q)^5 - 6(2-q)^5 & 1 \le q < 2 \\ (3-q)^5 & 2 \le q < 3 \\ 0 & q \ge 3\end{cases}\)	\(3.0\)	\(1/120\)	\(7/(478\pi)\)	\(1/(120\pi)\)
Wendland C2	'wendland_c2'	\(\begin{cases}(1-q/2)^4(2q+1) & q \le 2 \\ 0 & q > 2\end{cases}\)	\(2.0\)	\(5/8\)	\(7/(4\pi)\)	\(21/(16\pi)\)
Wendland C4	'wendland_c4'	\(\begin{cases}(1-q/2)^6\left(\frac{35}{12}q^2+3q+1\right) & q \le 2 \\ 0 & q > 2\end{cases}\)	\(2.0\)	\(3/4\)	\(9/(4\pi)\)	\(495/(256\pi)\)
Wendland C6	'wendland_c6'	\(\begin{cases}(1-q/2)^8(4q^3+6.25q^2+4q+1) & q \le 2 \\ 0 & q > 2\end{cases}\)	\(2.0\)	\(55/64\)	\(78/(28\pi)\)	\(1365/(512\pi)\)