@@ -60,29 +60,30 @@ struct MyGaussianKernel <: TrixiParticles.AbstractSmoothingKernel{2} end
6060
6161# By looking at the implementation of existing kernels in TrixiParticles.jl,
6262# we can see that a kernel implementation requires three functions.
63- # `TrixiParticles.kernel `, which is the kernel function itself,
64- # `TrixiParticles.kernel_deriv `, which is the derivative of the kernel function,
65- # and `TrixiParticles.compact_support`, which defines the compact support of the
66- # kernel in relation to the smoothing length.
63+ # `TrixiParticles.kernel_unsafe `, which is the kernel function itself,
64+ # `TrixiParticles.kernel_deriv_div_r_unsafe `, which is the derivative of the
65+ # kernel divided by ``r``, and `TrixiParticles.compact_support`, which defines
66+ # the compact support of the kernel in relation to the smoothing length.
6767# The latter is relevant for determining the search radius of the neighborhood search.
68- function TrixiParticles. kernel (kernel:: MyGaussianKernel , r, h)
68+ #
69+ # We implement `kernel_deriv_div_r_unsafe` instead of `kernel_deriv` directly since
70+ # this avoids an extra division in the hot loop and is robust near ``r=0``.
71+ # The public function `TrixiParticles.kernel_deriv` is defined automatically by
72+ # TrixiParticles from this method (and multiplies by ``r`` again when needed).
73+ # In the unsafe functions, we do not check the compact support; this is handled in the
74+ # safe wrappers `TrixiParticles.kernel` and `TrixiParticles.kernel_deriv` based on the
75+ # compact support defined in `TrixiParticles.compact_support`.
76+ function TrixiParticles. kernel_unsafe (kernel:: MyGaussianKernel , r:: Real , h)
6977 q = r / h
7078
71- if q < 2
72- return 1 / (pi * h^ 2 ) * exp (- q^ 2 )
73- end
74-
75- return 0.0
79+ return 1 / (pi * h^ 2 ) * exp (- q^ 2 )
7680end
7781
78- function TrixiParticles. kernel_deriv (kernel:: MyGaussianKernel , r, h)
82+ function TrixiParticles. kernel_deriv_div_r_unsafe (kernel:: MyGaussianKernel , r:: Real , h)
7983 q = r / h
8084
81- if q < 2
82- return 1 / (pi * h^ 2 ) * (- 2 * q) * exp (- q^ 2 ) / h
83- end
84-
85- return 0.0
85+ kernel_deriv = 1 / (pi * h^ 2 ) * (- 2 * q) * exp (- q^ 2 ) / h
86+ return kernel_deriv / r
8687end
8788
8889TrixiParticles. compact_support (:: MyGaussianKernel , h) = 2 * h
0 commit comments