Add all the other typically used smoothing kernels for SPH

src/general/smoothing_kernels.jl

@@ -6,11 +6,17 @@ of dimensions.
 
 We provide the following smoothing kernels:
 
-| Smoothing Kernel                          | Compact Support   |
-| :---------------------------------------- | :---------------- |
-| [`SchoenbergCubicSplineKernel`](@ref)     | $[0, 2h]$         |
-| [`SchoenbergQuarticSplineKernel`](@ref)   | $[0, 2.5h]$       |
-| [`SchoenbergQuinticSplineKernel`](@ref)   | $[0, 3h]$         |
+| Smoothing Kernel                          | Compact Support   | Typ. Smoothing Length | Recommended Application | Stability |
+| :---------------------------------------- | :---------------- | :-------------------- | :---------------------- | :-------- |
+| [`SchoenbergCubicSplineKernel`](@ref)     | $[0, 2h]$         | $1.1 to 1.3$          | General + sharp waves   | ++        |
+| [`SchoenbergQuarticSplineKernel`](@ref)   | $[0, 2.5h]$       | $1.1 to 1.5$          | General                 | +++       |
+| [`SchoenbergQuinticSplineKernel`](@ref)   | $[0, 3h]$         | $1.1 to 1.5$          | General                 | ++++      |
+| [`GaussianKernel`](@ref)                  | $[0, 3h]$         | $1.0 to 1.5$          | Literature              | +++++     |
+| [`WendlandC2Kernel`](@ref)                | $[0, 1h]$         | $2.5 to 4.0$          | **Try first**           | ++++      |
+| [`WendlandC4Kernel`](@ref)                | $[0, 1h]$         | $3.0 to 4.5$          | General                 | +++++     |
+| [`WendlandC6Kernel`](@ref)                | $[0, 1h]$         | $3.5 to 5.0$          | General                 | +++++     |
+| [`Poly6Kernel`](@ref)                     | $[0, 1h]$         | $1.5 to 2.5$          | Literature              | +         |
+| [`SpikyKernel`](@ref)                     | $[0, 1h]$         | $1.5 to 3.0$          | Sharp corners + waves   | +         |
 
 !!! note "Usage"
     The kernel can be called as
@@ -67,15 +73,19 @@ W(r, h) = \frac{\sigma_d}{h^d} e^{-r^2/h^2}
 
 where ``d`` is the number of dimensions and
 
-- `` \sigma_2 = \frac{1}{\pi} `` for 2D
-- `` \sigma_3 = \frac{1}{\pi^{3/2}} `` for 3D
+- `` \sigma_2 = \frac{1}{\pi} `` for 2D,
+- `` \sigma_3 = \frac{1}{\pi^{3/2}} `` for 3D.
 
 This kernel function has an infinite support, but in practice,
 it's often truncated at a certain multiple of ``h``, such as ``3h``.
 
-In this implementation, the kernel is truncated at `3h`,
+In this implementation, the kernel is truncated at ``3h``,
 so this kernel function has a compact support of ``[0, 3h]``.
 
+The smoothing length is typically in range `[1.0 : 1.5] * particle_spacing`.
+
+For general information and usage see [`SmoothingKernel`](@ref).
+
 Note:
 This truncation makes this Kernel not conservative,
 which is beneficial in regards to stability but makes it less accurate.
@@ -134,6 +144,8 @@ For an analytic formula for higher order Schoenberg kernels, see (Monaghan, 1985
 The largest disadvantage of Schoenberg Spline Kernel is the rather non-smooth first derivative,
 which can lead to increased noise compared to other kernel variants.
 
+The smoothing length is typically in range `[1.1 : 1.3] * particle_spacing`.
+
 For general information and usage see [`SmoothingKernel`](@ref).
 
 ## References:
@@ -213,6 +225,9 @@ For an analytic formula for higher order Schoenberg kernels, see (Monaghan, 1985
 
 The largest disadvantage of Schoenberg Spline Kernel are the rather non-smooth first derivative,
 which can lead to increased noise compared to other kernel variants.
+
+The smoothing length is typically in range `[1.1 : 1.5] * particle_spacing`.
+
 For general information and usage see [`SmoothingKernel`](@ref).
 
 ## References:
@@ -306,6 +321,9 @@ For an analytic formula for higher order Schoenberg kernels, see (Monaghan, 1985
 
 The largest disadvantage of Schoenberg Spline Kernel are the rather non-smooth first derivative,
 which can lead to increased noise compared to other kernel variants.
+
+The smoothing length is typically in range `[1.1 : 1.5] * particle_spacing`.
+
 For general information and usage see [`SmoothingKernel`](@ref).
 
 ## References:
@@ -376,16 +394,16 @@ abstract type WendlandKernel{NDIMS} <: SmoothingKernel{NDIMS} end
 Wendland C2 kernel (Wendland, 1995), a piecewise polynomial function designed to have compact support and to
 be twice continuously differentiable everywhere. Given by
 
-\[ W(r, h) = \frac{1}{h^d} w(r/h) \]
+`` W(r, h) = \frac{1}{h^d} w(r/h) ``
 
 with
 
-\[
+```math
 w(q) = \sigma \begin{cases}
     (1 - q)^4 (4q + 1)    & \text{if } 0 \leq q < 1, \\
     0                     & \text{if } q \geq 1,
 \end{cases}
-\]
+```
 
 where `` d `` is the number of dimensions and `` \sigma `` is a normalization factor dependent on the dimension.
 The normalization factor `` \sigma `` is `` 40/7\pi `` in two dimensions or `` 21/2\pi `` in three dimensions.
@@ -395,40 +413,47 @@ This kernel function has a compact support of `` [0, h] ``.
 For a detailed discussion on Wendland functions and their applications in SPH, see (Dehnen & Aly, 2012).
 The smoothness of these functions is also the largest disadvantage as they lose details at sharp corners.
 
+The smoothing length is typically in range `[2.5 : 4.0] * particle_spacing`.
+
+For general information and usage see [`SmoothingKernel`](@ref).
+
 ## References:
 - Walter Dehnen & Hassan Aly.
-"Improving convergence in smoothed particle hydrodynamics simulations without pairing instability".
-In: Monthly Notices of the Royal Astronomical Society 425.2 (2012), pages 1068-1082.
-[doi: 10.1111/j.1365-2966.2012.21439.x](https://doi.org/10.1111/j.1365-2966.2012.21439.x)
+  "Improving convergence in smoothed particle hydrodynamics simulations without pairing instability".
+  In: Monthly Notices of the Royal Astronomical Society 425.2 (2012), pages 1068-1082.
+  [doi: 10.1111/j.1365-2966.2012.21439.x](https://doi.org/10.1111/j.1365-2966.2012.21439.x)
 
 - Holger Wendland.
  "Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree."
- In: Advances in computational Mathematics 4 (1995): 389-396.
+ In: Advances in computational Mathematics 4 (1995), pages 389-396.
  [doi: 10.1007/BF02123482](https://doi.org/10.1007/BF02123482)
-
 """
 struct WendlandC2Kernel{NDIMS} <: WendlandKernel{NDIMS} end
 
 @fastpow @inline function kernel(kernel::WendlandC2Kernel, r::Real, h)
     q = r / h
 
+    result = (1 - q)^4 * (4q + 1)
+
     # Zero out result if q >= 1
-    result = ifelse(q < 1, normalization_factor(kernel, h) * (1 - q)^4 * (4q + 1), zero(q))
+    result = ifelse(q < 1, normalization_factor(kernel, h) * result, zero(q))
 
     return result
 end
 
-@fastpow @inline function kernel_deriv(kernel::WendlandC2Kernel, r::Real, h)
+@fastpow @muladd @inline function kernel_deriv(kernel::WendlandC2Kernel, r::Real, h)
     inner_deriv = 1 / h
     q = r * inner_deriv
 
     q1_3 = (1 - q)^3
     q1_4 = (1 - q)^4
 
+    result = -4 * q1_3 * (4q + 1)
+    result = result + q1_4 * 4
+
     # Zero out result if q >= 1
     result = ifelse(q < 1,
-                    normalization_factor(kernel, h) * (-4 * q1_3 * (4q + 1) + q1_4 * 4) *
-                    inner_deriv, zero(q))
+                    normalization_factor(kernel, h) * result * inner_deriv, zero(q))
 
     return result
 end
@@ -442,30 +467,34 @@ end
 Wendland C4 kernel, a piecewise polynomial function designed to have compact support and to
 be four times continuously differentiable everywhere. Given by
 
-\[ W(r, h) = \frac{1}{h^d} w(r/h) \]
+`` W(r, h) = \frac{1}{h^d} w(r/h) ``
 
 with
 
-\[
+```math
 w(q) = \sigma \begin{cases}
     (1 - q)^6 * (35q^2 / 3 + 6q + 1)   & \text{if } 0 \leq q < 1, \\
     0                                  & \text{if } q \geq 1,
 \end{cases}
-\]
+```
 
 where `` d `` is the number of dimensions and `` \sigma `` is a normalization factor dependent
-on the dimension. The exact value of `` \sigma `` needs to be determined based on the dimension to ensure proper normalization.
+on the dimension. The normalization factor `` \sigma `` is `` 9 / \pi `` in two dimensions or `` 495 / 32\pi `` in three dimensions.
 
 This kernel function has a compact support of `` [0, h] ``.
 
 For a detailed discussion on Wendland functions and their applications in SPH, see (Dehnen & Aly, 2012).
 The smoothness of these functions is also the largest disadvantage as they loose details at sharp corners.
 
+The smoothing length is typically in range `[3.0 : 4.5] * particle_spacing`.
+
+For general information and usage see [`SmoothingKernel`](@ref).
+
 ## References:
 - Walter Dehnen & Hassan Aly.
-"Improving convergence in smoothed particle hydrodynamics simulations without pairing instability".
-In: Monthly Notices of the Royal Astronomical Society 425.2 (2012), pages 1068-1082.
-[doi: 10.1111/j.1365-2966.2012.21439.x](https://doi.org/10.1111/j.1365-2966.2012.21439.x)
+  "Improving convergence in smoothed particle hydrodynamics simulations without pairing instability".
+  In: Monthly Notices of the Royal Astronomical Society 425.2 (2012), pages 1068-1082.
+  [doi: 10.1111/j.1365-2966.2012.21439.x](https://doi.org/10.1111/j.1365-2966.2012.21439.x)
 
 - Holger Wendland.
  "Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree."
@@ -477,10 +506,10 @@ struct WendlandC4Kernel{NDIMS} <: WendlandKernel{NDIMS} end
 @fastpow @inline function kernel(kernel::WendlandC4Kernel, r::Real, h)
     q = r / h
 
+    result = (1 - q)^6 * (35q^2 / 3 + 6q + 1)
+
     # Zero out result if q >= 1
-    result = ifelse(q < 1,
-                    normalization_factor(kernel, h) * (1 - q)^6 * (35q^2 / 3 + 6q + 1),
-                    zero(q))
+    result = ifelse(q < 1, normalization_factor(kernel, h) * result, zero(q))
 
     return result
 end
@@ -506,30 +535,34 @@ end
 
 Wendland C6 kernel, a piecewise polynomial function designed to have compact support and to be six times continuously differentiable everywhere. Given by:
 
-\[ W(r, h) = \frac{1}{h^d} w(r/h) \]
+`` W(r, h) = \frac{1}{h^d} w(r/h) ``
 
 with:
 
-\[
+```math
 w(q) = \sigma \begin{cases}
     (1 - q)^8 (32q^3 + 25q^2 + 8q + 1)    & \text{if } 0 \leq q < 1, \\
     0                                     & \text{if } q \geq 1,
 \end{cases}
-\]
+```
 
 where `` d `` is the number of dimensions and `` \sigma `` is a normalization factor dependent
-on the dimension. The exact value of `` \sigma `` needs to be determined based on the dimension to ensure proper normalization.
+on the dimension. The normalization factor `` \sigma `` is `` 78 / 7 \pi`` in two dimensions or `` 1365 / 64\pi`` in three dimensions.
 
 This kernel function has a compact support of `` [0, h] ``.
 
 For a detailed discussion on Wendland functions and their applications in SPH, see (Dehnen & Aly, 2012).
 The smoothness of these functions is also the largest disadvantage as they loose details at sharp corners.
 
+The smoothing length is typically in range `[3.5 : 5.0] * particle_spacing`.
+
+For general information and usage see [`SmoothingKernel`](@ref).
+
 ## References:
 - Walter Dehnen & Hassan Aly.
-"Improving convergence in smoothed particle hydrodynamics simulations without pairing instability".
-In: Monthly Notices of the Royal Astronomical Society 425.2 (2012), pages 1068-1082.
-[doi: 10.1111/j.1365-2966.2012.21439.x](https://doi.org/10.1111/j.1365-2966.2012.21439.x)
+  "Improving convergence in smoothed particle hydrodynamics simulations without pairing instability".
+  In: Monthly Notices of the Royal Astronomical Society 425.2 (2012), pages 1068-1082.
+  [doi: 10.1111/j.1365-2966.2012.21439.x](https://doi.org/10.1111/j.1365-2966.2012.21439.x)
 
 - Holger Wendland.
  "Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree."
@@ -541,10 +574,10 @@ struct WendlandC6Kernel{NDIMS} <: WendlandKernel{NDIMS} end
 @fastpow @inline function kernel(kernel::WendlandC6Kernel, r::Real, h)
     q = r / h
 
+    result = (1 - q)^8 * (32q^3 + 25q^2 + 8q + 1)
+
     # Zero out result if q >= 1
-    result = ifelse(q < 1,
-                    normalization_factor(kernel, h) * (1 - q)^8 * (32q^3 + 25q^2 + 8q + 1),
-                    zero(q))
+    result = ifelse(q < 1, normalization_factor(kernel, h) * result, zero(q))
 
     return result
 end
@@ -570,19 +603,19 @@ end
 
 Poly6 kernel, a commonly used kernel in SPH literature, especially in computer graphics contexts. It is defined as
 
-\[ W(r, h) = \frac{1}{h^d} w(r/h) \]
+`` W(r, h) = \frac{1}{h^d} w(r/h) ``
 
 with
 
-\[
+```math
 w(q) = \sigma \begin{cases}
     (1 - q^2)^3    & \text{if } 0 \leq q < 1, \\
     0              & \text{if } q \geq 1,
 \end{cases}
-\]
+```
 
 where `` d `` is the number of dimensions and `` \sigma `` is a normalization factor that depends
-on the dimension. The exact value of `` \sigma `` needs to be determined based on the dimension to ensure proper normalization.
+on the dimension. The normalization factor `` \sigma `` is `` 4 / \pi`` in two dimensions or `` 315 / 64\pi`` in three dimensions.
 
 This kernel function has a compact support of `` [0, h] ``.
 
@@ -591,6 +624,10 @@ other kernels that might offer better accuracy for hydrodynamic simulations. Fur
 its derivatives are not that smooth, which can lead to stability problems.
 It is also susceptible to clumping.
 
+The smoothing length is typically in range `[1.5 : 2.5] * particle_spacing`.
+
+For general information and usage see [`SmoothingKernel`](@ref).
+
 ## References:
 - Matthias Müller, David Charypar, and Markus Gross. "Particle-based fluid simulation for interactive applications".
   In: Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation. Eurographics Association. 2003, pages 154-159.
@@ -601,8 +638,10 @@ struct Poly6Kernel{NDIMS} <: SmoothingKernel{NDIMS} end
 @inline function kernel(kernel::Poly6Kernel, r::Real, h)
     q = r / h
 
+    result = (1 - q^2)^3
+
     # Zero out result if q >= 1
-    result = ifelse(q < 1, normalization_factor(kernel, h) * (1 - q^2)^3, zero(q))
+    result = ifelse(q < 1, normalization_factor(kernel, h) * result, zero(q))
 
     return result
 end
@@ -611,10 +650,11 @@ end
     inner_deriv = 1 / h
     q = r * inner_deriv
 
+    result = -6 * q * (1 - q^2)^2
+
     # Zero out result if q >= 1
     result = ifelse(q < 1,
-                    -6 * q * (1 - q^2)^2 * normalization_factor(kernel, h) * inner_deriv,
-                    zero(q))
+                    result * normalization_factor(kernel, h) * inner_deriv, zero(q))
     return result
 end
 
@@ -629,26 +669,30 @@ end
 The Spiky kernel is another frequently used kernel in SPH, especially due to its desirable
 properties in preserving features near boundaries in fluid simulations. It is defined as:
 
-\[ W(r, h) = \frac{1}{h^d} w(r/h) \]
+`` W(r, h) = \frac{1}{h^d} w(r/h) ``
 
 with:
 
-\[
+```math
 w(q) = \sigma \begin{cases}
     (1 - q)^3    & \text{if } 0 \leq q < 1, \\
     0            & \text{if } q \geq 1,
 \end{cases}
-\]
+```
 
-where `` d `` is the number of dimensions and `` \sigma `` is a normalization factor, which
-depends on the dimension and ensures the kernel integrates to 1 over its support.
+where `` d `` is the number of dimensions and the normalization factor `` \sigma `` is `` 10 / \pi``
+in two dimensions or `` 15 / \pi`` in three dimensions.
 
 This kernel function has a compact support of `` [0, h] ``.
 
 The Spiky kernel is particularly known for its sharp gradients, which can help to preserve
 sharp features in fluid simulations, especially near solid boundaries.
 These sharp gradients at the boundary are also the largest disadvantage as they can lead to instability.
 
+The smoothing length is typically in range `[1.5 : 3.0] * particle_spacing`.
+
+For general information and usage see [`SmoothingKernel`](@ref).
+
 ## References:
 - Matthias Müller, David Charypar, and Markus Gross. "Particle-based fluid simulation for interactive applications".
   In: Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation. Eurographics Association. 2003, pages 154-159.
@@ -659,8 +703,10 @@ struct SpikyKernel{NDIMS} <: SmoothingKernel{NDIMS} end
 @inline function kernel(kernel::SpikyKernel, r::Real, h)
     q = r / h
 
+    result = (1 - q)^3
+
     # Zero out result if q >= 1
-    result = ifelse(q < 1, normalization_factor(kernel, h) * (1 - q)^3, zero(q))
+    result = ifelse(q < 1, normalization_factor(kernel, h) * result, zero(q))
 
     return result
 end
@@ -669,9 +715,10 @@ end
     inner_deriv = 1 / h
     q = r * inner_deriv
 
+    result = -3 * (1 - q)^2
+
     # Zero out result if q >= 1
-    result = ifelse(q < 1, -3 * (1 - q)^2 * normalization_factor(kernel, h) * inner_deriv,
-                    zero(q))
+    result = ifelse(q < 1, result * normalization_factor(kernel, h) * inner_deriv, zero(q))
 
     return result
 end

test/general/smoothing_kernels.jl

@@ -13,8 +13,8 @@ using QuadGK
         return 4 * π * integral_3d_radial
     end
 
-    # All smoothing kernels should integrate to something close to 1
-    # Don't show all kernel tests in the final overview
+    # All smoothing kernels should integrate to something close to 1.
+    # Don't show all kernel tests in the final overview.
     @testset verbose=false "Integral" begin
         # Treat the truncated Gaussian kernel separately
         @testset "GaussianKernel" begin