Geometric Density Average Force SPH (GDSPH)

• GDSPH is a formulation of SPH that computes interparticle forces using the geometric mean of densities, effectively suppressing spurious surface tension artifacts.
• It is derived from Hamilton’s principle and variational methods, providing rigorous convergence and enhanced numerical stability across discontinuities and complex interfaces.
• Benchmark simulations in fluid mechanics, astrophysics, and multimaterial flows demonstrate GDSPH’s superior performance in accuracy, mixing, and computational efficiency.

Geometric Density Average Force Smoothed Particle Hydrodynamics (GDSPH) is an advanced formulation of the Smoothed Particle Hydrodynamics (SPH) method that modifies the computation of interparticle forces by employing a geometric mean of densities. This approach fundamentally improves the handling of density discontinuities, suppresses spurious surface tension artifacts in multiphase and interface-dominated flows, and enhances accuracy and robustness for simulations spanning fluid mechanics, astrophysics, and soft matter. Historical development traces to both variational derivations of continuum-to-particle systems and targeted efforts to alleviate known weaknesses in traditional SPH, most notably the density-weighted bias in force calculations.

1. Theoretical Foundations and Systematic Derivation

The derivation of GDSPH starts from Hamilton’s principle of least action for continuum mechanics. The continuum Lagrangian is expressed in terms of the velocity and an internal energy function that depends on the regularized density field, $\tilde{\rho}_t(\xi) = (W_h * \mu_t)(\xi)$, where $W_h$ is a symmetric mollifier (e.g., Wendland or Gaussian kernel) and $\mu_t$ is the measure-valued mass distribution (Evers et al., 2015). Discretization steps replace $\mu_0$ with a sum of Dirac deltas, yielding particle-based equations.

A critical aspect in GDSPH is the order in which regularization and action variation are applied. When the variation is taken before regularizing, the pressure gradient term becomes:

$$\mathbf{F}_i = -\frac{1}{\tilde{\rho}_i} \frac{d}{d\rho}\left[\tilde{\rho}_i^2 \bar{F}'(\tilde{\rho}_i)\right] \nabla \tilde{\rho}_i$$

This analytic evaluation of the pressure gradient omits explicit pairwise pressure terms, in contrast to traditional SPH, which uses symmetrized sums. The choice encapsulated by a switching parameter $\theta$ produces two branches: $\theta = 1$ for the classical (pairwise, momentum-conserving) SPH; $\theta = 0$ for the GDSPH variant (Evers et al., 2015). Rigorous convergence in Wasserstein distance has been established, with strong results for multiphase, shock, and drag-dominated flows.

2. GDSPH Discretization and Force Formulation

GDSPH replaces the conventional arithmetic density weighting in force expressions with a geometric mean. The canonical form of the pressure-gradient force for particles $i$ and $j$ is:

$$f_{p,i} = -\sum_j m_j \frac{P_i + P_j}{\rho_i \rho_j} \nabla_i W_{ij}$$

or as implemented in graphics and simulation:

$$f_{p,i} = -\sum_j m_j \frac{P_i + P_j}{G_{ij}^2} \nabla W_{ij}, \qquad G_{ij} = \sqrt{\rho_i \rho_j}$$

The symmetric operator ensures conservation and improved consistency with continuum mechanical equations (Wadsley et al., 2017, Koschier et al., 2020, Jeyaraj et al., 2020). In stiff or discontinuous density regimes, this geometric averaging suppresses spurious tensile instability and oscillations, especially across strong density gradients and interfaces.

A general SPH symmetric derivative in GDSPH adopts $\phi=1$ (versus $\phi=\rho$ in TSPH):

$$\frac{\nabla A}{\rho} \approx \sum_j \frac{m_j}{\rho_i \rho_j} (A_i + A_j) \nabla_i W_{ij}$$

This structure reduces zero-order errors $E_0$ typically present in conventional SPH (Wissing et al., 2019). Complementary volume correction factors $(2\rho_i)/(\rho_i + \rho_j)$ appear in stabilization and dissipation schemes (Zheng et al., 2023).

3. Robustness Across Discontinuities, Interfaces, and Boundaries

GDSPH has demonstrated significant advantages in treating density discontinuities where standard kernel-weighted SPH density estimation fails (Ruiz-Bonilla et al., 2022). In planetary impact simulations, uncorrected discontinuities introduce surface-tension-like artificial forces, suppressing mixing or causing particle ejection. GDSPH, and associated density correction schemes—such as imbalance-statistics-based blending—restore physically accurate pressures and densities near interfaces.

Boundary handling is further enhanced by geometric formulations of the Shepard renormalization factor, converting volume integrals truncated at free surfaces or walls into analytically tractable surface integrals. This method enables consistent force and density computations even with sparse particle support near boundaries while maintaining computational efficiency (Calderon-Sanchez et al., 14 Jan 2025).

4. Applications: Fluids, Magnetohydrodynamics, and Multimaterial Simulations

The GDSPH approach has found broad adoption in astrophysics (e.g., Gasoline2, Changa codes (Wadsley et al., 2017, Wissing et al., 2019)), computer graphics, and planetary science. In magnetohydrodynamics, GDSPH robustly handles MRI turbulence, jet launching, and magnetic tower formation in stratified accretion disks. Key attributes include:

• Removal of leading errors and surface tension artifacts (Kelvin–Helmholtz, multiphase instabilities).
• Accurate magnetic field evolution with constrained divergence cleaning (Wissing et al., 2019, Wissing et al., 2021).
• Efficient jet launching and magnetic structure formation under high dynamic range/resolution.
• Superior convergence rates for collapse and mixing-driven problems, compared to traditional SPH and mesh-based methods.

In compressible and incompressible fluids, GDSPH supports advanced solvers (predictive-corrective, divergence-free SPH), operator splitting, and strong coupling between fluids and solids (Koschier et al., 2020). The use of geometric density averaging accommodates:

• Pressure Poisson equations with improved stability.
• Enhanced material interface tracking in multimaterial flows.
• Large time steps in implicit schemes.

5. Computational and Algorithmic Considerations

Implementation of GDSPH introduces minimal computational overhead relative to TSPH. The only major change is the altered gradient formulation. When integrated with compact hashing neighborhood search, modern kernel (e.g., Wendland), gradient-based shock detection, and boundary corrections, GDSPH supports scalability to millions of interacting particles with robust accuracy (Koschier et al., 2020, Calderon-Sanchez et al., 14 Jan 2025).

Computational efficiency is maintained in three dimensions by analytical handling of divergent kernel terms near boundaries. The overall cost for boundary correction is typically $<10\%$ of total step time in high-resolution scenarios (Calderon-Sanchez et al., 14 Jan 2025). In multi-phase flows with large density ratios, volume correction factors regulate dissipation, suppress pressure oscillations, and stabilize particle distributions (Zheng et al., 2023). GDSPH also synergizes with generalized density dissipation and Riemann-solver-inspired stabilization techniques, allowing fine control of numerical diffusion.

6. Benchmark Results and Convergence Properties

GDSPH has been validated against a wide suite of benchmarks:

• Square and Kelvin–Helmholtz tests: Retains shape, supports instability growth, and minimizes spurious surface tension (Wadsley et al., 2017).
• Evrard collapse and magnetized cloud collapse: Accurate shock capture and jet formation; fast convergence at moderate resolution (Wissing et al., 2019).
• MRI simulations: Reproduces turbulence and dynamo action; avoids runaway field amplification; maintains low divergence error $\varepsilon_{\mathrm{divB}} \lesssim 10^{-2}$ (Wissing et al., 2021).
• Graphics and planetary impacts: Enhances mixing and interface reconstruction for high-resolution, multimaterial simulations (Koschier et al., 2020, Ruiz-Bonilla et al., 2022).
• Rayleigh–Taylor, bubble rising, and sloshing benchmarks: Suppresses pressure oscillations, ensures stability and accurate wave evolution, demonstrates agreement with experimental data (Zheng et al., 2023).

Mathematical convergence in the Wasserstein metric $\mathcal{W}(\mu_0^n, \mu_0) \to 0$, with uniform convergence in time for evolved measures (Evers et al., 2015), underpins the reliability of GDSPH.

7. Extensions, Limitations, and Future Directions

GDSPH is extensible to solid mechanics via geometric-average-based density/stress computations and corrected deformation gradients. Combining GDSPH with advanced boundary handling (Shepard factor geometric correction), stabilization and dissipation schemes offers improved robustness and accuracy for massive multiphysics simulations. A plausible implication is greater integration into commercial solvers and open-source codes as demands for multimaterial, interface-accurate simulation increase.

Limitations remain in density estimation at sharply discontinuous interfaces; recent density-correction approaches ameliorate these issues but may require case-specific tuning (Ruiz-Bonilla et al., 2022). For extreme density contrasts, geometric averaging must be carefully paired with volume corrections to prevent numerical artifacts.

Ongoing research trends point toward hybrid meshless-mesh approaches, adaptive kernel support scaling, and tighter coupling with Riemann-solver-based stabilization, all of which pair naturally with GDSPH’s foundational principles.

---

GDSPH represents a rigorously derived and comprehensively validated SPH variant that significantly improves accuracy and robustness in simulations involving discontinuities, multimaterial interfaces, and complex boundary geometries. Its geometric averaging of densities in force computations sets a new standard for momentum conservation, mixing, and stability across a range of computational physics domains.