LudwigBoess / SPHKernels.jl

Julia implementation of some common SPH kernels
MIT License
10 stars 0 forks source link

Simplify kernel structs, value and derivative for different dimensions. #14

Closed LudwigBoess closed 2 years ago

LudwigBoess commented 3 years ago

Also as discussed in #11 the kernel structs need to be reworked.

As a first step we'll drop the n_neighbours entry, as it is not very useful. Second the dimensionality can be stored in the kernel to unify a lot of functions. For example for the Cubic kernel:

struct Cubic{T} <: SPHKernel
    dim::Int64
    norm::T
end

function Cubic(T::DataType=Float64, dim::Integer=3)

    if dim == 1
        norm = 4.0/3.0
    elseif dim == 2 
        norm = 40.0/7π 
    elseif dim == 3
        norm = 8.0/π
    else
        error("Cubic not defined for $dim dimensions!")
    end

    return Cubic{T}(dim, norm)

end

"""
    kernel_value(kernel::Cubic{T}, u::Real, h_inv::Real) where T

Evaluate cubic spline at position ``u = \\frac{x}{h}``.
"""
function kernel_value(kernel::Cubic{T}, u::Real, h_inv::Real) where T

    n = kernel.norm * h_inv^kernel.dim

    if u < 0.5
        return ( 1 + 6 * ( u - 1 ) * u^2) * n |> T
    elseif u < 1
        u_m1 = 1 - u 
        return  2u_m1^3 * n |> T
    else
        return 0.0 |> T
    end

end

This only raises the question how to handle kernels that have different functional forms for different dimensions, e.g. the Wendland kernels.

One option would be to implement a proper type hirarchy, something like:

abstract type AbstractSPHKernel end

abstract type WendlandC6 <: AbstractSPHKernel end

and then have different constructors and let multiple dispatch take care of the rest.

AhmedSalih3d commented 3 years ago

Hi!

Thanks for splitting it into sub-issues, we can have much more detailed discussions this way.

I am not really good at this part of structs etc., but I can see your idea here. You want to defined each individual kernel to be part of an upper-class "AbstractSPHKernel". By doing it this way, we also associate the dim with the kernel it self, which means that it only has to be defined once and we let go of W1,W2 principle (which the user them selves can redefine if useful for them I suppose).

Could you explain this point to me though:

This only raises the question how to handle kernels that have different functional forms for different dimensions, e.g. the Wendland kernels.

The Wendland kernel I have been using has one form for all dimensions (only defined for 2 and 3) but the constant norm changes as for any other kernel. I suppose you are talking about a kernel you use a lot your self which has the same name or?

My Concluding Remarks

I think this is clearly the right direction. Removing any extra values (n_neighbours) is very smart and makes it very clear what the purpose of the structs are. If you have a clear idea of where to take this, I would say go for it, since it will no matter what make the outcome better.

Kind regards

LudwigBoess commented 2 years ago

Hi Ahmed,

please have a look at #17. Do you think that makes sense? I really like the change, it's easier to maintain and more difficult to introduce bugs... Turns out there was actually another typo in WendlandC8 that I only just found after I simplified it...

Cheers, Ludwig

AhmedSalih3d commented 2 years ago

Hi again

Just had a look at it, it looks really nice! I wondered if it is possible for me to download the pr code you made in some way, I just can't find it in Github right now?

So I had to read through the code and it looked to me that you did exactly as we agreed on, so I am looking forward to trying it.

And no problem about the typo's, don't worry too much! Everyone using public libraries has their own responsibility for verifying what they use, works, so don't feel like everything is on your "shoulders" :) Atleast we (you) ironed out one more mistake!

Kind regards

LudwigBoess commented 2 years ago

You should be able to clone the repo and then use git checkout simplify_kernel_functions to switch to the development branch.

LudwigBoess commented 2 years ago

Since we agreed that this is the way forward I merged this now into the roadmap_to_2_0 branch. So you need to use that one if you want to play around with it! :)

I think that resolves this issue, so I'll close it now.