Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Simplify kernel functions #17

Merged
merged 5 commits into from
Sep 11, 2021
Merged

Simplify kernel functions #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0adfdd8
simplify kernels (possible fix for #14)
LudwigBoess Sep 10, 2021
bb4268d
updated docs for new kernel definition
LudwigBoess Sep 10, 2021
e2bc9d7
renamed SPHKernel to AbstractSPHKernel for more consistent naming
LudwigBoess Sep 10, 2021
da99eae
tests for error handling of undefined dimensions
LudwigBoess Sep 10, 2021
9ff66df
removed sph functions (handled in #18)
LudwigBoess Sep 10, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 59 additions & 11 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@

This package supplies a number of kernels frequently used in Smoothed-Particle Hydrodynamics (SPH), as well as functions to evaluate their values and derivatives in 2D and 3D.

These kernels include the B-splines (`Cubic` and `Quintic`) suggested in [Monaghan & Lattanzio (1985)](https://ui.adsabs.harvard.edu/abs/1985A%26A...149..135M/abstract) and the Wendland functions (`WendlandC2`, `WendlandC4` and `WendlandC6`) as suggested in [Dehnen & Aly (2012)](https://academic.oup.com/mnras/article/425/2/1068/1187211).
These kernels include the B-splines (`Cubic` and `Quintic`) suggested in [Monaghan & Lattanzio (1985)](https://ui.adsabs.harvard.edu/abs/1985A%26A...149..135M/abstract) and the Wendland functions (`WendlandC2`, `WendlandC4`, `WendlandC6` and `WendlandC8` from [Wendland (2009)](https://www.researchgate.net/publication/220179293_Divergence-Free_Kernel_Methods_for_Approximating_the_Stokes_Problem)) as suggested in [Dehnen & Aly (2012)](https://academic.oup.com/mnras/article/425/2/1068/1187211).

In this implementation we follow the convention of Dehnen&Aly in using the 'compact kernel support' as a means to define the maximum extent of the kernel. They denote this ``H`` in their paper, for convenience (aka for not having to type caps) we use `h` in the code.

Expand All @@ -15,29 +15,77 @@ In this implementation we follow the convention of Dehnen&Aly in using the 'comp
To evaluate a 3D kernel you need to use the function

```julia
kernel_value_3D(k::SPHKernel, u::Real, h_inv::Real)
kernel_value(k::AbstractSPHKernel, u::Real, h_inv::Real)
```

where `SPHKernel` is the supertype for an implemented SPH kernel, ``u = \frac{x}{h}`` is the distance to the kernel origin in measures of the compact kernel support and `h_inv` is the inverse of the compact kernel support.
where [AbstractSPHKernel](@ref) is the supertype for an implemented SPH kernel, ``u = \frac{x}{h}`` is the distance to the kernel origin in measures of the compact kernel support and `h_inv` is the inverse of the compact kernel support.

The same goes for a 1D kernel

```julia
kernel_value_1D(k::SPHKernel, u::Real, h_inv::Real)
kernel_value(k::AbstractSPHKernel, u::Real, h_inv::Real)
```

and a 2D kernel

If you want your code to look a little more fancy you can also use the alternative functions [𝒲₁](@ref), where the respective subscript refers to the dimension:

```julia
kernel_value_2D(k::SPHKernel, u::Real, h_inv::Real)
𝒲( kernel::AbstractSPHKernel, u::Real, h_inv::Real) = kernel_value(kernel, u, h_inv)
```

If you want your code to look a little more fancy you can also use the alternative functions `𝒲₁`, where the respective subscript refers to the dimension:
As an example:
```julia
using SPHKernels

# Wendland C6 kernel with double precision in 3D
k = WendlandC6(Float64, 3)
# distance between the particle and the origin of the kernel
r = 0.5
h = 1.0
h_inv = 1.0/h
u = r * h_inv

# kernel value at position r
val = 𝒲(k, u, h_inv)
```


## Evaluating Derivatives

Similar to [Evaluating Kernels](@ref) you can evluate a kernel derivative with

```julia
kernel_deriv(k::AbstractSPHKernel, u::Real, h_inv::Real)
```

or in the fancy way:

```julia
𝒲₁( kernel::SPHKernel, u::Real, h_inv::Real ) = kernel_value_1D(kernel, u, h_inv)
𝒲₂( kernel::SPHKernel, u::Real, h_inv::Real ) = kernel_value_2D(kernel, u, h_inv)
𝒲₃( kernel::SPHKernel, u::Real, h_inv::Real ) = kernel_value_3D(kernel, u, h_inv)
d𝒲(kernel::AbstractSPHKernel, u::Real, h_inv::Real) = kernel_deriv(kernel, u, h_inv)
```

Please see the docs for more details!
## Bias Correction

You can correct for the kernel bias of the Wendland kernels as described in [Dehnen & Aly (2012)](https://academic.oup.com/mnras/article/425/2/1068/1187211), Eq. 18 + 19 with the functions:

```julia
bias_correction(kernel::AbstractSPHKernel, density::Real, m::Real, h_inv::Real, n_neighbours::Integer)
```

or again in the fancy way

```julia
δρ(kernel::AbstractSPHKernel, density::Real, m::Real, h_inv::Real, n_neighbours::Integer) = bias_correction(kernel, density, m, h_inv, n_neighbours)

```

This will return a new value for the density:

```julia
using SPHKernels
density = 1.0
kernel = WendlandC6(3)

# correct density
density = bias_correction_3D(kernel, density, 1.0, 0.5, 295)
```
30 changes: 18 additions & 12 deletions docs/src/extending.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,32 +7,38 @@ DocTestSetup = quote
end
```

If you need a different kernel function than the ones I implemented you can add them by defining a new kernel `struct` as a subtype of [SPHKernel](@ref)
If you need a different kernel function than the ones I implemented you can add them by defining a new kernel `struct` as a subtype of [AbstractSPHKernel](@ref)

```@example 1
using SPHKernels # hide
struct MyKernel{T} <: SPHKernel
n_neighbours::Int64
norm_1D::T
norm_2D::T
norm_3D::T
struct MyKernel{T} <: AbstractSPHKernel
dim::Int64
norm::T
end

"""
MyKernel(T::DataType=Float64, n_neighbours::Integer=295)
MyKernel(T::DataType=Float64, dim::Integer=3)

Set up a `MyKernel` kernel for a given DataType `T`.
"""
MyKernel(T::DataType=Float64, n_neighbours::Integer=1000) = MyKernel{T}(n_neighbours, T(1), T(1), T(1))
MyKernel(T::DataType=Float64, dim::Integer=3) = MyKernel{T}(dim, T(1))


"""
MyKernel(dim::Integer)

Define `MyKernel` kernel with dimension `dim` for the native `DataType` of the OS.
"""
MyKernel(dim::Integer) = MyKernel{T}(typeof(1.0), dim)
```

and defining its value and derivative in 2D and 3D, e.g.
and defining its value and derivative, e.g.

```@example 1
function kernel_value_2D(kernel::MyKernel{T}, u::Real, h_inv::Real) where T
function kernel_value(kernel::MyKernel{T}, u::Real, h_inv::Real) where T

if u < 1
n = kernel.norm_2D * h_inv^2
n = kernel.norm * h_inv^kernel.dim
return 1.0 * n |> T
else
return 0.0 |> T
Expand All @@ -42,7 +48,7 @@ end
```

```@example 1
@inline function kernel_deriv_2D(kernel::MyKernel{T}, u::Real, h_inv::Real) where T
@inline function kernel_deriv(kernel::MyKernel{T}, u::Real, h_inv::Real) where T
return 0.0 |> T
end
```
Expand Down
42 changes: 15 additions & 27 deletions docs/src/kernels.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ end

This package supplies a number of kernels frequently used in Smoothed-Particle Hydrodynamics (SPH), as well as functions to evaluate their values and derivatives in 2D and 3D.

These kernels include the B-splines ([Cubic](@ref) and [Quintic](@ref)) suggested in [Monaghan & Lattanzio (1985)](https://ui.adsabs.harvard.edu/abs/1985A%26A...149..135M/abstract) and the Wendland functions ([WendlandC2](@ref), [WendlandC4](@ref) and [WendlandC6](@ref)) as suggested in [Dehnen & Aly (2012)](https://academic.oup.com/mnras/article/425/2/1068/1187211).
These kernels include the B-splines ([Cubic](@ref) and [Quintic](@ref)) suggested in [Monaghan & Lattanzio (1985)](https://ui.adsabs.harvard.edu/abs/1985A%26A...149..135M/abstract) and the Wendland functions ([WendlandC2](@ref), [WendlandC4](@ref), [WendlandC6](@ref)) and [WendlandC8](@ref) ([Wendland 2009](https://www.researchgate.net/publication/220179293_Divergence-Free_Kernel_Methods_for_Approximating_the_Stokes_Problem)) as suggested in [Dehnen & Aly (2012)](https://academic.oup.com/mnras/article/425/2/1068/1187211).

In this implementation we follow the convention of Dehnen&Aly in using the 'compact kernel support' as a means to define the maximum extent of the kernel. They denote this ``H`` in their paper, for convenience (aka for not having to type caps) we use `h` in the code.

Expand All @@ -18,43 +18,38 @@ In this implementation we follow the convention of Dehnen&Aly in using the 'comp
To evaluate a 3D kernel you need to use the function

```julia
kernel_value_3D(k::SPHKernel, u::Real, h_inv::Real)
kernel_value(k::AbstractSPHKernel, u::Real, h_inv::Real)
```

where [SPHKernel](@ref) is the supertype for an implemented SPH kernel, ``u = \frac{x}{h}`` is the distance to the kernel origin in measures of the compact kernel support and `h_inv` is the inverse of the compact kernel support.
where [AbstractSPHKernel](@ref) is the supertype for an implemented SPH kernel, ``u = \frac{x}{h}`` is the distance to the kernel origin in measures of the compact kernel support and `h_inv` is the inverse of the compact kernel support.

The same goes for a 1D kernel

```julia
kernel_value_1D(k::SPHKernel, u::Real, h_inv::Real)
kernel_value(k::AbstractSPHKernel, u::Real, h_inv::Real)
```

and a 2D kernel

```julia
kernel_value_2D(k::SPHKernel, u::Real, h_inv::Real)
```

If you want your code to look a little more fancy you can also use the alternative functions [𝒲₁](@ref), where the respective subscript refers to the dimension:

```julia
𝒲₁( kernel::SPHKernel, u::Real, h_inv::Real) = kernel_value_1D(kernel, u, h_inv)
𝒲₂( kernel::SPHKernel, u::Real, h_inv::Real) = kernel_value_2D(kernel, u, h_inv)
𝒲₃( kernel::SPHKernel, u::Real, h_inv::Real) = kernel_value_3D(kernel, u, h_inv)
𝒲( kernel::AbstractSPHKernel, u::Real, h_inv::Real) = kernel_value(kernel, u, h_inv)
```

As an example:
```@example
using SPHKernels # hide
k = WendlandC6()

# Wendland C6 kernel with double precision in 3D
k = WendlandC6(Float64, 3)
# distance between the particle and the origin of the kernel
r = 0.5
h = 1.0
h_inv = 1.0/h
u = r * h_inv

# kernel value at position r
val = 𝒲₃(k, u, h_inv)
val = 𝒲(k, u, h_inv)

println("val = $val")
```
Expand All @@ -65,34 +60,27 @@ println("val = $val")
Similar to [Evaluating Kernels](@ref) you can evluate a kernel derivative with

```julia
kernel_deriv_3D(k::SPHKernel, u::Real, h_inv::Real)
kernel_deriv(k::AbstractSPHKernel, u::Real, h_inv::Real)
```

or in the fancy way:

```julia
d𝒲₁(kernel::SPHKernel, u::Real, h_inv::Real) = kernel_deriv_1D(kernel, u, h_inv)
d𝒲₂(kernel::SPHKernel, u::Real, h_inv::Real) = kernel_deriv_2D(kernel, u, h_inv)
d𝒲₃(kernel::SPHKernel, u::Real, h_inv::Real) = kernel_deriv_3D(kernel, u, h_inv)

d𝒲(kernel::AbstractSPHKernel, u::Real, h_inv::Real) = kernel_deriv(kernel, u, h_inv)
```

## Bias Correction

You can correct for the kernel bias of the Wendland kernels as described in [Dehnen & Aly (2012)](https://academic.oup.com/mnras/article/425/2/1068/1187211), Eq. 18 + 19 with the functions:

```julia
bias_correction_1D(kernel::SPHKernel, density::Real, m::Real, h_inv::Real)
bias_correction_2D(kernel::SPHKernel, density::Real, m::Real, h_inv::Real)
bias_correction_3D(kernel::SPHKernel, density::Real, m::Real, h_inv::Real)
bias_correction(kernel::AbstractSPHKernel, density::Real, m::Real, h_inv::Real, n_neighbours::Integer)
```

or again in the fancy way

```julia
δρ₁(kernel::SPHKernel, density::Real, m::Real, h_inv::Real) = bias_correction_1D(kernel, density, m, h_inv)
δρ₂(kernel::SPHKernel, density::Real, m::Real, h_inv::Real) = bias_correction_2D(kernel, density, m, h_inv)
δρ₃(kernel::SPHKernel, density::Real, m::Real, h_inv::Real) = bias_correction_3D(kernel, density, m, h_inv)
δρ(kernel::AbstractSPHKernel, density::Real, m::Real, h_inv::Real, n_neighbours::Integer) = bias_correction(kernel, density, m, h_inv, n_neighbours)

```

Expand All @@ -101,10 +89,10 @@ This will return a new value for the density:
```@example
using SPHKernels # hide
density = 1.0
kernel = WendlandC6()
kernel = WendlandC6(3)

# correct density
density = bias_correction_3D(kernel, density, 1.0, 0.5)
density = bias_correction_3D(kernel, density, 1.0, 0.5, 295)

println("density = $density")
```
76 changes: 12 additions & 64 deletions src/SPHKernels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,19 +4,10 @@

module SPHKernels

export kernel_value_1D, 𝒲₁,
kernel_value_2D, 𝒲₂,
kernel_value_3D, 𝒲₃,
kernel_deriv_1D, d𝒲₁,
kernel_deriv_2D, d𝒲₂,
kernel_deriv_3D, d𝒲₃,
kernel_gradient_1D, ∇𝒲₁,
kernel_gradient_2D, ∇𝒲₂,
kernel_gradient_3D, ∇𝒲₃,
bias_correction_1D, δρ₁,
bias_correction_2D, δρ₂,
bias_correction_3D, δρ₃,
SPHKernel,
export kernel_value, 𝒲,
kernel_deriv, d𝒲,
bias_correction, δρ,
AbstractSPHKernel,
Cubic,
Quintic,
WendlandC2,
Expand All @@ -26,85 +17,42 @@ module SPHKernels


"""
SPHKernel
AbstractSPHKernel

Supertype for all SPH kernels.
"""
abstract type SPHKernel end
abstract type AbstractSPHKernel end

include("bsplines/Cubic.jl")
include("bsplines/Quintic.jl")
include("wendland/C2.jl")
include("wendland/C4.jl")
include("wendland/C6.jl")
include("wendland/C8.jl")
include("sph_functions/gradients.jl")

# multiple dispatch for nicer look

"""
𝒲₁( kernel::SPHKernel, u::Real, h_inv::Real)
𝒲₁( kernel::AbstractSPHKernel, u::Real, h_inv::Real)

Evaluate 1D spline at position ``u = \\frac{x}{h}``.
"""
𝒲₁( kernel::SPHKernel, u::Real, h_inv::Real) = kernel_value_1D(kernel, u, h_inv)

"""
d𝒲₁( kernel::SPHKernel, u::Real, h_inv::Real)

Evaluate 1D derivative at position ``u = \\frac{x}{h}``.
"""
d𝒲₁(kernel::SPHKernel, u::Real, h_inv::Real) = kernel_deriv_1D(kernel, u, h_inv)

"""
𝒲₂( kernel::SPHKernel, u::Real, h_inv::Real)

Evaluate 2D spline at position ``u = \\frac{x}{h}``.
"""
𝒲₂( kernel::SPHKernel, u::Real, h_inv::Real) = kernel_value_2D(kernel, u, h_inv)

"""
d𝒲₂( kernel::SPHKernel, u::Real, h_inv::Real)

Evaluate 1D derivative at position ``u = \\frac{x}{h}``.
"""
d𝒲₂(kernel::SPHKernel, u::Real, h_inv::Real) = kernel_deriv_2D(kernel, u, h_inv)
𝒲( kernel::AbstractSPHKernel, u::Real, h_inv::Real) = kernel_value(kernel, u, h_inv)

"""
𝒲₃( kernel::SPHKernel, u::Real, h_inv::Real)

Evaluate 3D spline at position ``u = \\frac{x}{h}``.
"""
𝒲₃( kernel::SPHKernel, u::Real, h_inv::Real) = kernel_value_3D(kernel, u, h_inv)

"""
d𝒲₃( kernel::SPHKernel, u::Real, h_inv::Real)
d𝒲₂( kernel::AbstractSPHKernel, u::Real, h_inv::Real)

Evaluate 1D derivative at position ``u = \\frac{x}{h}``.
"""
d𝒲₃(kernel::SPHKernel, u::Real, h_inv::Real) = kernel_deriv_3D(kernel, u, h_inv)
d𝒲(kernel::AbstractSPHKernel, u::Real, h_inv::Real) = kernel_deriv(kernel, u, h_inv)

"""
δρ₁(kernel::SPHKernel, density::Real, m::Real, h_inv::Real)
δρ₁(kernel::AbstractSPHKernel, density::Real, m::Real, h_inv::Real)

Corrects the 1D density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
"""
δρ₁(kernel::SPHKernel, density::Real, m::Real, h_inv::Real) = bias_correction_1D(kernel, density, m, h_inv)
δρ(kernel::AbstractSPHKernel, density::Real, m::Real, h_inv::Real, n_neighbours::Integer) = bias_correction(kernel, density, m, h_inv, n_neighbours)

"""
δρ₂(kernel::SPHKernel, density::Real, m::Real, h_inv::Real)

Corrects the 2D density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
"""
δρ₂(kernel::SPHKernel, density::Real, m::Real, h_inv::Real) = bias_correction_2D(kernel, density, m, h_inv)

"""
δρ₃(kernel::SPHKernel, density::Real, m::Real, h_inv::Real)

Corrects the 3D density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
"""
δρ₃(kernel::SPHKernel, density::Real, m::Real, h_inv::Real) = bias_correction_3D(kernel, density, m, h_inv)



end # module
Loading