Linear-Corrected Gradient SPH (ISPH)

• Linear-corrected gradient SPH (ISPH) is a meshless method that applies local linear corrections to accurately recover gradients of linear functions.
• It enhances kernel normalization and partition of unity, reducing interpolation errors crucial for capturing fluid instabilities and complex flows.
• ISPH maintains conservation of mass, momentum, and energy through a variational formulation, ensuring reliable performance in multi-physics simulations.

Linear-Corrected Gradient SPH (ISPH) denotes a family of smoothed particle hydrodynamics (SPH) methods in which the approximation of differential operators—especially spatial gradients—is systematically improved via local linear correction. These methods are formulated to achieve improved consistency (the exact recovery of gradients of linear functions) while maintaining essential conservation properties, a capability of critical importance for robust fluid dynamics and astrophysical simulations.

1. Mathematical Foundations and Integral Gradient Formulation

Linear-corrected gradient SPH relies on replacing the SPH gradient approximation with an integral representation. For a scalar function $f(\mathbf{r})$, the gradient at position $\mathbf{r}$ is written as

$$\mathbf{I}(\mathbf{r}) = \int_{V} [f(\mathbf{r}') - f(\mathbf{r})] (\mathbf{r}' - \mathbf{r}) W(|\mathbf{r}' - \mathbf{r}|, h)\, \mathrm{d}^3\mathbf{r}'\,.$$

Expanding $f(\mathbf{r}')$ via a first-order Taylor expansion gives

$$f(\mathbf{r}') - f(\mathbf{r}) = \nabla f \cdot (\mathbf{r}' - \mathbf{r}) + \mathcal{O}(|\mathbf{r}' - \mathbf{r}|^2)\,,$$

leading to a matrix relation

$[\partial_i f] = \tau^{-1}_{ij} I_j\,,$

where

$$\tau_{ij} = \int (\mathbf{r}'_i - \mathbf{r}_i)(\mathbf{r}'_j - \mathbf{r}_j) W(|\mathbf{r}' - \mathbf{r}|, h)\, \mathrm{d}^3\mathbf{r}'\,,$$

and

$$I_j = \int [f(\mathbf{r}') - f(\mathbf{r})] (\mathbf{r}'_j - \mathbf{r}_j) W(|\mathbf{r}' - \mathbf{r}|, h) \, \mathrm{d}^3\mathbf{r}'\,.$$

In the discrete SPH setting, these integrals are evaluated as sums over neighboring particles, yielding element-wise forms per particle (indexed by $a$):

$$\tau_{ij,a} = \sum_{b} \frac{m_b}{\rho_b} (x_{i,b} - x_{i,a})(x_{j,b} - x_{j,a}) W_{ab}(h_a)$$

$$I_{j,a} = \sum_{b} m_b (x_{j,b} - x_{j,a}) W_{ab}(h_a)$$

The linear-corrected gradient estimate at particle $a$ is then:

$$\left( \frac{\partial f}{\partial x_i} \right)_a = \sum_{j=1}^d c_{ij,a} I_{j,a}$$

with $C = \tau^{-1}$. This construction is often termed IAD or IAD$_0$ in the literature (Garcia-Senz et al., 2011, Valdarnini, 2016, García-Senz et al., 2021).

2. Conservative Formulation and Variational Structure

A distinguishing property of the ISPH family is that, despite adopting a corrected gradient, the fundamental equations of motion and thermodynamics can be consistently derived from a discrete Lagrangian, ensuring exact conservation of mass, linear momentum, angular momentum, and energy (Garcia-Senz et al., 2011, Hopkins, 2012, García-Senz et al., 2021). The equations of motion assume forms such as

$$\frac{d^2 x_{i,a}}{dt^2} = - \sum_{b} m_b \left[ \frac{P_a}{\Omega_a \rho_a^2} \mathcal{A}_{i,ab}(h_a) + \frac{P_b}{\Omega_b \rho_b^2} \mathcal{A}'_{i,ab}(h_b) \right]$$

with

$$\mathcal{A}_{i,ab}(h_a) = \sum_{j=1}^d c_{ij,a}(x_{j,b} - x_{j,a}) W_{ab}(h_a)$$

and analogous expressions for $\mathcal{A}'$ at particle $b$. The factor $\Omega$ is related to the grad-$h$ correction and established from the density’s dependence on the smoothing length.

This variational derivation is essential: it guarantees compatibility with the conservative SPH infrastructure and ensures that modifications introduced for accuracy do not violate fundamental physical constraints (Garcia-Senz et al., 2011, García-Senz et al., 2021).

3. Gradient Correction, Partition of Unity, and Volume Elements

ISPH’s core accuracy benefit is contingent on both the local kernel normalization and the partition of unity. The volume element $V_a$ in the SPH summations is pivotal; traditionally, $V_a = m_a/\rho_a$, but generalized volume elements of the form

$$V_a = \frac{X_a}{k_a}, \quad k_a = \sum_b X_b W_{ab}(h_a)$$

with $X_a = m_a / \rho_a^0$ and $\rho_a^0$ as the standard SPH density estimate, are shown to greatly improve kernel normalization and partition of unity:

• The normalization error $E_1 = | \sum_b V_b W_{ab} - 1 |$ is minimized.
• The linear consistency error $E_2 = | \sum_b V_b (r_b - r_a) W_{ab} | / h_a$ is reduced.

The improved partition of unity realized through these volume elements directly enhances gradient estimation, particularly in problems with strong density contrasts (García-Senz et al., 2021).

4. Comparison with Standard and Alternative SPH Methods

Extensive benchmark studies reveal that linear-corrected (integral approach) SPH exhibits significantly lower errors and better physical fidelity in challenging regimes:

• For linear density profiles, the relative gradient error is much smaller than standard SPH.
• Growth of Kelvin–Helmholtz and Rayleigh–Taylor instabilities is more robust and less susceptible to noise or suppression, even with subsonic initial perturbations (Garcia-Senz et al., 2011, Valdarnini, 2016).
• In turbulence and vortex-driven problems, the inertial range and spectral properties reproduce theoretical expectations more faithfully than standard SPH (Valdarnini, 2016, Rosswog, 2014).
• Compared to advanced meshless and moving-mesh schemes, ISPH yields comparable accuracy and error norms, sometimes outperforming in the subsonic and high-shear limits (Valdarnini, 2016).

A summary comparison:

Gradient Type	Conservation	Linear Consistency	Robustness (subsonic/shear)
Standard (kernel)	Yes	No	Low
ISPH (integral)	Yes	Approximate/Yes	High
Renormalization (KGC)	Depends	Yes	Case-dependent

ISPH achieves improved interpolation accuracy with only moderate computational overhead (local $d \times d$ matrix inversion per particle, $d=2,3$).

5. Extensions to Physical Applications and Testing

ISPH methodology has been validated in various physically relevant test cases:

• Kelvin–Helmholtz and Rayleigh–Taylor Instabilities: ISPH triggers the correct instability growth even with minimal initial perturbations, whereas standard SPH fails (Garcia-Senz et al., 2011, Valdarnini, 2016).
• Sun-like Polytrope: 3D equilibrium simulations of self-gravitating polytropes show improved equilibrium structures and more accurate balance between pressure and gravity, with better conservation of the center-of-mass motion—although relaxation to equilibrium can be slightly slower than SPH (Garcia-Senz et al., 2011).
• Kelvin–Helmholtz, Cloud–Wind, and Sedov Explosion: ISPH accurately preserves interface shapes and instability dynamics, with better entropy evolution and reduced unphysical “blobs” compared to standard density-based SPH (García-Senz et al., 2021).
• Keplerian Disc and Turbulence: Angular momentum transport is better controlled, and turbulence power spectra more closely follow theoretical predictions, which is essential for astrophysical and galactic dynamics simulations (Valdarnini, 2016).

6. Connections to Modern Consistency Corrections and Theoretical Implications

Recent developments, such as the reverse kernel gradient correction (RKGC) (Zhang et al., 28 May 2024, Zhang et al., 14 Mar 2025), underscore the theoretical trajectory initiated by ISPH for reconciling high-order consistency with exact conservation. RKGC achieves exact zero- and first-order consistency in an anti-symmetric, conservative framework. While ISPH integral corrections are Lagrangian-conservative and can be tailored for high accuracy, methods such as RKGC suggest pathways to even higher-order consistency, e.g., via “reverse” KGC application, albeit with the challenge of robust particle relaxation to fully enforce the consistency conditions.

This line of research illustrates the broad impact of linear- (or higher-order) corrected gradients in ensuring that SPH is competitive with mesh-based approaches for both fluid and structure problems. The ISPH approach and its descendants remain a central theme in SPH algorithmic accuracy and robustness for engineering, astrophysics, and multi-physics FSI contexts, particularly under large density gradients or strong shear (García-Senz et al., 2021, Zhang et al., 28 May 2024, Zhang et al., 14 Mar 2025).

7. Current Challenges, Limitations, and Prospective Developments

While ISPH and successors provide clear accuracy and conservation advantages, certain challenges remain:

• Particle Disorder and Relaxation: High-order consistency corrections often assume (or are optimal for) well-relaxed particle distributions, which require robust pre-processing or dynamic relaxation algorithms.
• Interface Instability and Tensile Effects: Handling sharp density contrasts may necessitate surgical deviations from strict Lagrangian formulations (e.g., via a local $\sigma$ parameter) to suppress the tensile instability at interfaces, balancing conservation and robustness (García-Senz et al., 2021).
• Kernel Choice: The adoption of higher-order kernels (e.g., Wendland) can further improve gradient robustness and suppress noise but may require larger neighbor numbers for comparable accuracy (Rosswog, 2014, Valdarnini, 2016).
• Free Surface and Boundary Effects: Corrected gradient methods require careful handling at boundaries and free surfaces, as standard correction matrices can become ill-conditioned in highly non-uniform regions (Zhang et al., 14 Mar 2025).
• Extension to Full Incompressible Flows: While ISPH is dominant in both compressible and incompressible regimes (e.g., via pressure-Poisson or projection methods), maintaining stability and unique solvability in the incompressible limit hinges on particle connectivity, regularity, and time step constraints (Imoto, 2018, Muta et al., 2019).

Ongoing research targets more robust volume element definitions, improved dynamic stability in multi-physics settings, and further extension of high-order anti-symmetric consistent operators for both fluids and solid mechanics.

---

Linear-corrected gradient SPH (ISPH), anchored in the integral approach, has become a foundational methodology in meshless Lagrangian hydrodynamics modeling. Through direct enhancement of operator consistency and normalization, combined with careful compatibility with conservation laws, ISPH establishes a reliable, physically faithful, and extensible framework for simulation of complex, multi-physics flows.