LudwigBoess
commented
3 years ago Hi Ahmed,
thank you for your nice message and your good work on implementing a different kernel!
To adress your questions and points:
- grad(W): Yes, you are correct! Calling
∇𝒲the derivative is sloppy and inconsistent. I usually compute the actual gradient seperately, but it would make sense to supply an additional function to do it, as you did. Should we rename∇𝒲₃tod𝒲₃(and the same for the other dimentsions) ? - n_neighbours: I originally implemented this as a reference point. In the Gadget code you can define a range number to define how far the actual numbers of neighbours can deviate from the defined number. I have not made use of that for my work though, so it could also be dropped.
- Disabling the if statement: I think that would make it not very readable, so I'd avoid that. Unless you have a good suggestion how to do that elegantly? As you already noted the if statement uses up barely any of the compute time.
Please note that I defined the kernels so far to have u = r/h where r is the euclidean distance and h is the compact kernel support, as you already mentioned in your other comment. This should be in the range x ∈ [0, 1]. Your kernel seems to use x ∈ [0, 2].
Other then that please feel free to implement the kernels you want in your forked repo and then issue a pull request to the main repo! :)
I found some performance improvements for your gradient computation:
function kernel_grad_2D(kernel::DSPH_Quintic{T}, u::Real, h_inv::Real, x::AbstractArray, y::AbstractArray, r::Real) where T
n = kernel.norm_2D * h_inv^2
if u < 2
u_d1 = 5/8
u_d2 = (u - 2)
dWdq = u_d1 * u * (u_d2 * u_d2 * u_d2) |> T
else
return 0.0, 0.0 .|> T
end
# Calculate result for each direction
x_diff = x[1] - y[1]
y_diff = x[2] - y[2]
tmp = dWdq*r*h_inv
Wgx = tmp*x_diff |> T
Wgy = tmp*y_diff |> T
# returning as Tuple instead of Array avoids 1 allocation
# doing the type conversion only once gets rid of the other allocation
return Wgx, Wgy
end
# BenchmarkTools:
@btime Wgx, Wgy = kernel_grad_2D($k, $u2, $h_inv,$x,$y,$r)
# 4.870 ns (0 allocations: 0 bytes) I like your idea of having an actual gradient implementation in the package! Maybe it would be more convenient to have a more general implementation, like so:
function kernel_grad_2D(k::SPHKernel, h_inv::Real, pos1::Vector{<:Real}, pos2::Vector{<:Real})
x_diff = pos1[1] - pos2[1]
y_diff = pos1[2] - pos2[2]
r = √(x_diff^2 + y_diff^2)
u = r * h_inv
dW = kernel_deriv_2D(k, u, h_inv)
tmp = dW * u
Wgx = tmp * x_diff
Wgy = tmp * y_diff
return Wgx, Wgy
endWhat do you think?
AhmedSalih3d
So we should have a function, which takes any arbitrary AbstractSPHKernel and allows to calculate these properties, for example the gradient based on the dW given in the struct. Syntax perhaps like:
Hi!
Thank you very much for this effort. I have spend quite a bit of time on this code to produce some nice examples of what I want to ask about, so I hope you enjoy the read :)
Basically I had to import a kernel I know just to be sure about what I am doing - I have looked at how you code and tried to replicate it as close as possible - seems to also improve performance which is very nice, benchmarks at the bottom. Extending the code was quite easy to, you have made a very nice setup for this kind of thing.
I don't agree with your notation saying that grad(W) = kernel_deriv_dim and you can see the full explanation of why in kernel_grad_2d. What are your thoughts about this?
If you agree I think it should be split in kernel_value, kernel_deriv (dWdq) and kernel_gradient. Else I am looking forward to learn if I have misunderstood something; which could also be likely.
The code should run instantly without any other dependencies than the library here.
Kind regards