Add gradient, divergence and rotation functions

[LudwigBoess]

Currently ∇𝒲 is the same as kernel_derivative, which is of course wrong.

In a future release this will be fixed to d𝒲.

In the same spirit it would be useful to implement functions for gradient, divergence and rotation. How should we go about this?

One option Ahmed suggested for gradient computation is:

"""
    kernel_grad_3D(k::SPHKernel, h_inv::Real, part_pos::Vector{<:Real}, neighbour_pos::Vector{<:Real})

Computes the 3D gradient of the kernel `k` at the position of the neighbour `neighbour_pos`. 

``∇𝒲(x_i) = \\frac{d𝒲}{dx}\\vert_{x_j} \\frac{Δx_{ij}}{||x_{ij}||} \\frac{1}{h_i}``
"""
function kernel_grad_3D(k::SPHKernel, h_inv::Real, part_pos::Vector{<:Real}, neighbour_pos::Vector{<:Real})

    x_diff = part_pos[1] - neighbour_pos[1]
    y_diff = part_pos[2] - neighbour_pos[2]
    z_diff = part_pos[3] - neighbour_pos[3]

    r = √(x_diff^2 + y_diff^2 + z_diff^2)
    u = r * h_inv

    tmp = d𝒲₃(k, u, h_inv) / r
    Wgx = tmp*x_diff
    Wgy = tmp*y_diff
    Wgz = tmp*z_diff

    return Wgx, Wgy, Wgz
end

# Benchmark:
# 12.866 ns (0 allocations: 0 bytes)

This returns a Tuple which makes it quite performant. The downside is that it's dimension dependent, something we should move away from (see #14).

Defining it on an Array basis would help, but is less performant:

"""
    kernel_grad(k::SPHKernel, h_inv::Real, part_pos::Vector{<:Real}, neighbour_pos::Vector{<:Real})

Computes the gradient of the kernel `k` at the position of the neighbour `neighbour_pos`. 

``∇𝒲(x_i) = \\frac{d𝒲}{dx}\\vert_{x_j} \\frac{Δx_{ij}}{||x_{ij}||} \\frac{1}{h_i}``
"""
function kernel_grad(k::SPHKernel, h_inv::Real, part_pos::Vector{<:Real}, neighbour_pos::Vector{<:Real})

    Δx = part_pos - neighbour_pos

    r = √( sum( Δx.^2 ) )
    u = r * h_inv

    dW_r = d𝒲₃(k, u, h_inv) / r

    return Δx .* dW_r
end

# Benchmark
# 85.307 ns (3 allocations: 336 bytes)

Any suggestions?

[LudwigBoess]

Same goes for rotations, except that it's only defined in 3D.

To to compute the rotation of a quantity Q we could use something like:

function kernel_rotation(k::SPHKernel, h_inv::Real, 
                         part_pos::Vector{<:Real}, neighbour_pos::Vector{<:Real},
                         part_Q::Vector{<:Real}, neighbour_Q::Vector{<:Real})

    Δx = part_pos[1] - neighbour_pos[1]
    Δy = part_pos[2] - neighbour_pos[2]
    Δz = part_pos[3] - neighbour_pos[3]

    ΔQx = part_Q[1] - neighbour_Q[1]
    ΔQy = part_Q[2] - neighbour_Q[2]
    ΔQz = part_Q[3] - neighbour_Q[3]

    r = √(Δx^2 + Δy^2 + Δz^2)
    u = r * h_inv

    dWk_r = d𝒲₃(k, u, h_inv) / r

    return  dWk_r * ( Δy * ΔQz - Δz * ΔQy), 
            dWk_r * ( Δz * ΔQx - Δx * ΔQz), 
            dWk_r * ( Δx * ΔQy - Δy * ΔQx)

end

# Benchmark
# 15.042 ns (0 allocations: 0 bytes)

Or alternatively:

using LinearAlgebra

function kernel_rotation(k::SPHKernel, h_inv::Real, 
                         part_pos::Vector{<:Real}, neighbour_pos::Vector{<:Real},
                         part_Q::Vector{<:Real}, neighbour_Q::Vector{<:Real})

    Δx = part_pos - neighbour_pos
    ΔQ = part_Q - neighbour_Q

    r = √( sum( Δx.^2 ) )
    u = r * h_inv

    dW_r = d𝒲₃(k, u, h_inv) / r

    return  dWk_r .* (  Δx \times ΔQ )

end

# Benchmark
# 152.898 ns (5 allocations: 560 bytes)

[AhmedSalih3d]

Interesting work and thanks for iniatiating it!

Fixing the name of kernel_rotation In general maths and CFD my experience has been that this quantity is coined "curl". Would you mind too much that it would be named this? This means we would have:

kernel_value kernel_derivative kernel_gradient kernel_curl

Which I think is much clearer, since it follows the naming of common math operators.

Tuple's versus Arrays I think we should really considering getting it to work with the allocation less approach. It will make a difference when at scale and since Julia is often known for fast performance, I see us two (and hopefully more in the future!) as those who have to struggle a bit more to get it to work, to make it easier for others :-)

One idea I had, which I will mention now (but cannot fully code it due to time constraint today) is that perhaps Tullio.jl, could sort this problem out for us. It basically allows to use Einstein notation in Julia, which means that perhaps even if the input is an array, we can return it as a tuple. The size of the point vector input would then be implicitly known, since @tullio p[i], should loop through it. Perhaps "eachindex" could be used as well.

In my experience tuple has been the fastest to use and better than StaticArrays (if only performance) is considered. From an aesthetic view I also like the first function more in kernel_rotation, I like using "smallest operators" possible, but I prefer performance more - in this case both of these are fulfilled.

Sorry for not being able to contribute that much code as of now, but really nice to see this work here.

Kind regards

[LudwigBoess]

I opened #18 to handle development

[AhmedSalih3d]

@LudwigBoess to avoid dimension dependendness give this approach a shot:

function kernel_gradient(kernel_derivative,h_inv,r,p1,p2)
    q    = r*h_inv
    dWdq = kernel_derivative(q)

    tup = map(p1,p2) do i,j
        d = i - j
        dWdq*(d/r)*h_inv
    end

    return tup
end

I asked on the discourse and was advised to do it like this - I was playing with ForwardDiff again. Here I did the OOP you do not like with kernel_derivative, so you can skip that part, the key part is in the "tup = map". Basically it takes p1 and p2 as a tuple (should work for vector too) of any shape and returns the equivivalent gradient term. The full code for my example is shown below:

using BenchmarkTools
using ForwardDiff
using FastPow

const h      = 0.028284271247
const h_inv  = 1/h
const aD    = (7.0)/(4.0*π*h^2)

@fastpow Wendland(q) = aD*(1-q/2)^4 * (2*q+1)
@fastpow WendlandDerivative(q) = aD*((5/8)*q*(q-2)^3)

const df    = x -> ForwardDiff.derivative(Wendland, x)

@btime WendlandDerivative($(Ref(1.1))[]);
@btime df($(Ref(1.1))[]);

const pts = rand(1000);

@btime WendlandDerivative.(u) setup=(u=pts);

@btime df.(u) setup=(u=pts);

function kernel_gradient(kernel_derivative,h_inv,r,p1,p2)
    q    = r*h_inv
    dWdq = kernel_derivative(q)

    tup = map(p1,p2) do i,j
        d = i - j
        dWdq*(d/r)*h_inv
    end

    return tup
end

# Test Stuff Out

p1 = (1.0,1.0)
p2 = (1.0,1.025)
r  = (p1 .- p2).^2 |> sum |> sqrt

@btime kernel_gradient($WendlandDerivative,$h_inv,$r,$p1,$p2)
@btime kernel_gradient($df,$h_inv,$r,$p1,$p2)

# Implementation is too "clever", have to use "$Ref(p1)[]" etc. or more
# elements to compare the two.

# I generate these vectors in a bad way, but you see my point I hope.. 
p1a = Vector{Tuple}()
p2a   = Vector{Tuple}()
for _ in 1:1000
    tmp = tuple(rand(),rand(),rand())
    push!(p2a,tmp)
end

for _ in 1:1000
    push!(p1a,p2a[1])
end

ra  = Vector{Float64}()
for i in 1:1000
    push!(ra,(p1a[1] .- p2a[i]).^2 |> sum |> sqrt)
end
ra .+= 1e-6 #Avoid NaN-Values, small numerical error.

@btime kernel_gradient.($WendlandDerivative,$h_inv,$ra,$p1a,$p2a);
@btime kernel_gradient.($df,$h_inv,$ra,$p1a,$p2a);

What I like about this is that I can specify any kernel_derivative function and it will evaluate it ÀND that if it is 1D, 2D or 3D points, will not matter - this means I do not have to rely on multiple dispatch for kernel_gradient(kernel::Type), but only have one function I need to maintain. The ForwardDiff aspect is mostly for fun, the results for the last two benchmarks were:

Just wanted to share, since the map approach might be of interest to your implementation.

Kind regards

[LudwigBoess]

Hi Ahmed,

thanks for the code example with tup = map that looks like a very good solution! :)

Concerning the ForwardDiff, I can't reproduce your results. Taking the first 15 lines of your example and deleting the const (since in an actual run these won't be constants) I get:

@btime WendlandDerivative($(Ref(1.1))[])
#    1.919 ns (0 allocations: 0 bytes)

@btime df($(Ref(1.1))[])
#   20.848 ns (2 allocations: 32 bytes)

compared to the current implementation:

@btime d𝒲₃($(Ref(k))[], $(Ref(0.5))[], $(Ref(2.0))[])
#   3.006 ns (0 allocations: 0 bytes)

The fact that the current implementation is slower mainly comes from your version doing fewer operations since you already have the norm defined, I guess.

Also please consider that the numeric differentiation is just an approximation, so you loose precision:

WendlandDerivative(1.1) == df(1.1)
# gives false in my test

As impressive as I find the ForwardDiff package, I don't think it's actually useful here. The whole purpose of these kernels are that they approximate a Gaussian as best possible and have analytic, simple solutions for at least the first 2 derivatives. Having these derivatives defeats the purpose of having and autodiff packge do them for you, in my opinion.

[AhmedSalih3d]

Hi again, no problem at all regarding map, hopefully it works :) Let me know if there is any issue with that point.

ForwardDiff was just extra and not something we need to use - my main point was the map functionality.

The reason I put const infront of different stuff, was that I was evaluating in global scope. If everything was inside a function the performance of both approaches should be nearly equal.

Yep, it might not give the correct result to 1e-6 etc. When doing floating points comparison I would resort to isapprox and not strict equality too, which then should have returned true.

So don't worry, I agree with you that ForwardDiff is not something to be added directly to this package, since using it once would be a bit lazy hehe :-)

Kind regards

[LudwigBoess]

Hi Ahmed,

I finally had some time to look into this. Could you please check #22 and tell me what you think?

Cheers, Ludwig

[AhmedSalih3d]

Hi Ludwig!

Hope you are well.

I gave it a very short look now and can have a bit deeper look on Sunday if you want me to. In general it looks to me that you have managed to implement the concept in a great way, similar to what we were discussing i.e. having a unified interface to take the gradient etc. independent of the selected kernel. I also really like your use of unicode, even though "x_j" would be basically the same as: "xⱼ", I believe that using this since we can in Julia makes the code much easier readable and maintable. Also easier to debug functions this way in my experience (just a side note).

I noted in the docs that some LaTeX equations were not being rendered properly: https://github.com/LudwigBoess/SPHKernels.jl/blob/roadmap_to_2_0/docs/src/sph.md I believe I have seen a fix for this somewhere, I will let you know when I find it.

I am a bit new to test suite in Julia, but I suppose it would make sense for me to try to install the roadmap_to_2 version of the repos and play around with the new code from a user perspective and check validity of results too. This might be in a week or two though, I hope that is fine.

It is really starting to shape up to something very promising.

Kind regards, Ahmed

[LudwigBoess]

Hi Ahmed,

thank you for your kind feedback. Yes, if you have the time it would be great if you could check how it works from the user end (a bit of a sanity check if you will).

The LaTeX should look fine once Documenter.jl renders to .md files to html. See the rendered docs if the current version here.

Cheers, Ludwig

[LudwigBoess]

Hi Ahmed,

I decided to merge the branch now, as I need the updated dependency for some of my other packages. There's a disclaimer now on the version number problem, so people should sanity-check the results for now. I did my best in testing, but maybe I missed something.

I'm closing this issue for now.

Cheers, Ludwig

[AhmedSalih3d]

Hi Ahmed,

I decided to merge the branch now, as I need the updated dependency for some of my other packages. There's a disclaimer now on the version number problem, so people should sanity-check the results for now. I did my best in testing, but maybe I missed something.

I'm closing this issue for now.

Cheers, Ludwig

Hi again!

No worries, sorry for forgetting to get back to you. I got/am really busy, so it left my mind. Happy to see that you are using it in your other work.

Kind regards, Ahmed