diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 62acdcbba..031bccca3 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -13,7 +13,7 @@ using OrdinaryDiffEq fluid_particle_spacing = 0.02 # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model -boundary_layers = 3 +boundary_layers = 4 spacing_ratio = 1 boundary_particle_spacing = fluid_particle_spacing / spacing_ratio @@ -39,8 +39,8 @@ tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fl # ========================================================================================== # ==== Fluid -smoothing_length = 1.2 * fluid_particle_spacing -smoothing_kernel = SchoenbergCubicSplineKernel{2}() +smoothing_length = 3.0 * fluid_particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() fluid_density_calculator = ContinuityDensity() viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 1b1c88c5e..44be7fc42 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -41,7 +41,8 @@ export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, - SchoenbergQuinticSplineKernel + SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel, + WendlandC6Kernel, SpikyKernel, Poly6Kernel export StateEquationIdealGas, StateEquationCole export ArtificialViscosityMonaghan, ViscosityAdami export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, diff --git a/src/general/smoothing_kernels.jl b/src/general/smoothing_kernels.jl index 2abb0e29d..32f353bc1 100644 --- a/src/general/smoothing_kernels.jl +++ b/src/general/smoothing_kernels.jl @@ -6,11 +6,20 @@ 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$ | General (recommended) | ++++ | +| [`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 | + | + +We recommend to use the [`WendlandC2Kernel`](@ref) for most applications. +If less smoothing is needed, try [`SchoenbergCubicSplineKernel`](@ref), for more smoothing try [`WendlandC6Kernel`](@ref). !!! note "Usage" The kernel can be called as @@ -57,6 +66,62 @@ end kernel_correction_coefficient(system, particle) end +@doc raw""" + GaussianKernel{NDIMS}() + +Gaussian kernel given by +```math +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. + +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``, +so this kernel function has a compact support of ``[0, 3h]``. + +The smoothing length is typically in the range ``[1.0\delta, 1.5\delta]``, +where ``\delta`` is the typical 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. +""" +struct GaussianKernel{NDIMS} <: SmoothingKernel{NDIMS} end + +@inline @fastmath function kernel(kernel::GaussianKernel, r::Real, h) + q = r / h + + # Zero out result if q >= 3 + result = ifelse(q < 3, normalization_factor(kernel, h) * exp(-q^2), zero(q)) + + return result +end + +@inline @fastmath function kernel_deriv(kernel::GaussianKernel, r::Real, h) + inner_deriv = 1 / h + q = r * inner_deriv + + # Zero out result if q >= 3 + result = ifelse(q < 3, + -2 * q * normalization_factor(kernel, h) * exp(-q^2) * inner_deriv, + zero(q)) + + return result +end + +@inline compact_support(::GaussianKernel, h) = 3 * h + +@inline normalization_factor(::GaussianKernel{2}, h) = 1 / (pi * h^2) +@inline normalization_factor(::GaussianKernel{3}, h) = 1 / (pi^(3 / 2) * h^3) + @doc raw""" SchoenbergCubicSplineKernel{NDIMS}() @@ -80,6 +145,12 @@ This kernel function has a compact support of ``[0, 2h]``. For an overview of Schoenberg cubic, quartic and quintic spline kernels including normalization factors, see (Price, 2012). 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 the range ``[1.1\delta, 1.3\delta]``, +where ``\delta`` is the typical particle spacing. + For general information and usage see [`SmoothingKernel`](@ref). ## References: @@ -156,6 +227,13 @@ This kernel function has a compact support of ``[0, 2.5h]``. For an overview of Schoenberg cubic, quartic and quintic spline kernels including normalization factors, see (Price, 2012). 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 the range ``[1.1\delta, 1.5\delta]``, +where ``\delta`` is the typical particle spacing. + For general information and usage see [`SmoothingKernel`](@ref). ## References: @@ -246,6 +324,13 @@ This kernel function has a compact support of ``[0, 3h]``. For an overview of Schoenberg cubic, quartic and quintic spline kernels including normalization factors, see (Price, 2012). 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 the range ``[1.1\delta, 1.5\delta]``, +where ``\delta`` is the typical particle spacing. + For general information and usage see [`SmoothingKernel`](@ref). ## References: @@ -304,3 +389,364 @@ end @inline normalization_factor(::SchoenbergQuinticSplineKernel{1}, h) = 1 / 120h @inline normalization_factor(::SchoenbergQuinticSplineKernel{2}, h) = 7 / (478 * pi * h^2) @inline normalization_factor(::SchoenbergQuinticSplineKernel{3}, h) = 1 / (120 * pi * h^3) + +abstract type WendlandKernel{NDIMS} <: SmoothingKernel{NDIMS} end + +# Compact support for all Wendland kernels +@inline compact_support(::WendlandKernel, h) = h + +@doc raw""" + WendlandC2Kernel{NDIMS}() + +Wendland C2 kernel (Wendland, 1995), a piecewise polynomial function designed to have compact support and to +be twice continuously differentiable everywhere. Given by + +```math + 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. + +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 the range ``[2.5\delta, 4.0\delta]``, +where ``\delta`` is the typical 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) + +- Holger Wendland. + "Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree." + 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) * result, zero(q)) + + return result +end + +@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) * result * inner_deriv, zero(q)) + + return result +end + +@inline normalization_factor(::WendlandC2Kernel{2}, h) = 7 / (pi * h^2) +@inline normalization_factor(::WendlandC2Kernel{3}, h) = 21 / (2pi * h^3) + +@doc raw""" + WendlandC4Kernel{NDIMS}() + +Wendland C4 kernel, a piecewise polynomial function designed to have compact support and to +be four times continuously differentiable everywhere. Given by + +```math + 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 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 the range ``[3.0\delta, 4.5\delta]``, +where ``\delta`` is the typical 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) + +- Holger Wendland. + "Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree." + In: Advances in computational Mathematics 4 (1995): 389-396. + [doi: 10.1007/BF02123482](https://doi.org/10.1007/BF02123482) +""" +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) * result, zero(q)) + + return result +end + +@fastpow @muladd @inline function kernel_deriv(kernel::WendlandC4Kernel, r::Real, h) + q = r / h + term1 = (1 - q)^6 * (6 + 70 / 3 * q) + term2 = 6 * (1 - q)^5 * (1 + 6q + 35 / 3 * q^2) + derivative = term1 - term2 + + # Zero out result if q >= 1 + result = ifelse(q < 1, normalization_factor(kernel, h) * derivative / h, + zero(derivative)) + + return result +end + +@inline normalization_factor(::WendlandC4Kernel{2}, h) = 9 / (pi * h^2) +@inline normalization_factor(::WendlandC4Kernel{3}, h) = 495 / (32pi * h^3) + +@doc raw""" + WendlandC6Kernel{NDIMS}() + +Wendland C6 kernel, a piecewise polynomial function designed to have compact support and to be six times continuously differentiable everywhere. Given by: + +```math +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 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 the range ``[3.5\delta, 5.0\delta]``, +where ``\delta`` is the typical 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) + +- Holger Wendland. + "Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree." + In: Advances in computational Mathematics 4 (1995): 389-396. + [doi: 10.1007/BF02123482](https://doi.org/10.1007/BF02123482) +""" +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) * result, zero(q)) + + return result +end + +@fastpow @muladd @inline function kernel_deriv(kernel::WendlandC6Kernel, r::Real, h) + q = r / h + term1 = -8 * (1 - q)^7 * (32q^3 + 25q^2 + 8q + 1) + term2 = (1 - q)^8 * (96q^2 + 50q + 8) + derivative = term1 + term2 + + # Zero out result if q >= 1 + result = ifelse(q < 1, normalization_factor(kernel, h) * derivative / h, + zero(derivative)) + + return result +end + +@inline normalization_factor(::WendlandC6Kernel{2}, h) = 78 / (7pi * h^2) +@inline normalization_factor(::WendlandC6Kernel{3}, h) = 1365 / (64pi * h^3) + +@doc raw""" + Poly6Kernel{NDIMS}() + +Poly6 kernel, a commonly used kernel in SPH literature, especially in computer graphics contexts. It is defined as + +```math +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 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] ``. + +Poly6 is well-known for its computational simplicity, though it's worth noting that there are +other kernels that might offer better accuracy for hydrodynamic simulations. Furthermore, +its derivatives are not that smooth, which can lead to stability problems. +It is also susceptible to clumping. + +The smoothing length is typically in the range ``[1.5\delta, 2.5\delta]``, +where ``\delta`` is the typical 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. + [doi: 10.5555/846276.846298](https://doi.org/10.5555/846276.846298) +""" +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) * result, zero(q)) + + return result +end + +@inline function kernel_deriv(kernel::Poly6Kernel, r::Real, h) + 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, + result * normalization_factor(kernel, h) * inner_deriv, zero(q)) + return result +end + +@inline compact_support(::Poly6Kernel, h) = h + +@inline normalization_factor(::Poly6Kernel{2}, h) = 4 / (pi * h^2) +@inline normalization_factor(::Poly6Kernel{3}, h) = 315 / (64pi * h^3) + +@doc raw""" + SpikyKernel{NDIMS}() + +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: + +```math + 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 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 the range ``[1.5\delta, 3.0\delta]``, +where ``\delta`` is the typical 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. + [doi: 10.5555/846276.846298](https://doi.org/10.5555/846276.846298) +""" +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) * result, zero(q)) + + return result +end + +@inline function kernel_deriv(kernel::SpikyKernel, r::Real, h) + inner_deriv = 1 / h + q = r * inner_deriv + + result = -3 * (1 - q)^2 + + # Zero out result if q >= 1 + result = ifelse(q < 1, result * normalization_factor(kernel, h) * inner_deriv, zero(q)) + + return result +end + +@inline compact_support(::SpikyKernel, h) = h + +@inline normalization_factor(::SpikyKernel{2}, h) = 10 / (pi * h^2) +@inline normalization_factor(::SpikyKernel{3}, h) = 15 / (pi * h^3) diff --git a/test/Project.toml b/test/Project.toml index c18f9cde1..1ecb52bbf 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,6 +3,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] diff --git a/test/general/general.jl b/test/general/general.jl index 3947021f0..04d2dd5d1 100644 --- a/test/general/general.jl +++ b/test/general/general.jl @@ -1,5 +1,6 @@ # `util.jl` is added here to avoid confusion with `test_util.jl` include("util.jl") include("initial_condition.jl") +include("smoothing_kernels.jl") include("density_calculator.jl") include("semidiscretization.jl") diff --git a/test/general/smoothing_kernels.jl b/test/general/smoothing_kernels.jl new file mode 100644 index 000000000..87498f6d6 --- /dev/null +++ b/test/general/smoothing_kernels.jl @@ -0,0 +1,50 @@ +using QuadGK + +@testset verbose=true "Smoothing Kernels" begin + function integrate_kernel_2d(smk) + integral_2d_radial, _ = quadgk(r -> r * TrixiParticles.kernel(smk, r, 1.0), 0, + TrixiParticles.compact_support(smk, 1.0), rtol=1e-15) + return 2 * π * integral_2d_radial + end + + function integrate_kernel_3d(smk) + integral_3d_radial, _ = quadgk(r -> r^2 * TrixiParticles.kernel(smk, r, 1.0), 0, + TrixiParticles.compact_support(smk, 1.0), rtol=1e-15) + 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. + @testset verbose=false "Integral" begin + # Treat the truncated Gaussian kernel separately + @testset "GaussianKernel" begin + error_2d = abs(integrate_kernel_2d(GaussianKernel{2}()) - 1.0) + error_3d = abs(integrate_kernel_3d(GaussianKernel{3}()) - 1.0) + + # Large error due to truncation + @test 1e-4 < error_2d < 1e-3 + @test 1e-4 < error_3d < 1e-3 + end + + # All other kernels + kernels = [ + SchoenbergCubicSplineKernel, + SchoenbergQuarticSplineKernel, + SchoenbergQuinticSplineKernel, + WendlandC2Kernel, + WendlandC4Kernel, + WendlandC6Kernel, + SpikyKernel, + Poly6Kernel, + ] + + @testset "$kernel" for kernel in kernels + # The integral should be 1 for all kernels + error_2d = abs(integrate_kernel_2d(kernel{2}()) - 1.0) + error_3d = abs(integrate_kernel_3d(kernel{3}()) - 1.0) + + @test error_2d <= 1e-15 + @test error_3d <= 1e-15 + end + end +end