improve performance of kernel(::SchoenbergQuarticSplineKernel, r, h),…

… take 2 (#255)

Co-authored-by: Erik Faulhaber <[email protected]>

Project.toml

@@ -6,9 +6,11 @@ version = "0.1.0"
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
+FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Morton = "2a6d852e-3fac-5a38-885c-fe708af2d09e"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -21,8 +23,10 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 DiffEqCallbacks = "2.25"
+FastPow = "0.1"
 ForwardDiff = "0.10"
 Morton = "0.1"
+MuladdMacro = "0.2"
 Polyester = "0.7.5"
 Reexport = "1"
 SciMLBase = "1, 2"

src/TrixiParticles.jl

@@ -4,8 +4,10 @@ using Reexport: @reexport
 
 using Dates
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
+using FastPow: @fastpow
 using LinearAlgebra: norm, dot, I, tr
 using Morton: cartesian2morton
+using MuladdMacro: @muladd
 using Polyester: Polyester, @batch
 using Printf: @printf, @sprintf
 using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!,

src/general/smoothing_kernels.jl

@@ -163,24 +163,29 @@ For an analytic formula for higher order kernels, see (Monaghan, 1985).
 """
 struct SchoenbergQuarticSplineKernel{NDIMS} <: SmoothingKernel{NDIMS} end
 
-function kernel(kernel::SchoenbergQuarticSplineKernel, r::Real, h)
+# Note that `floating_point_number^integer_literal` is lowered to `Base.literal_pow`.
+# Currently, specializations reducing this to simple multiplications exist only up
+# to a power of three, see
+# https://github.com/JuliaLang/julia/blob/34934736fa4dcb30697ac1b23d11d5ad394d6a4d/base/intfuncs.jl#L327-L339
+# Higher powers just call a generic power implementation. Here, we accept to lose
+# some precision but gain performance by using plain multiplications instead via
+# `@fastpow`.
+@muladd @fastpow @inline function kernel(kernel::SchoenbergQuarticSplineKernel, r::Real, h)
     q = r / h
 
-    if q >= 5 / 2
-        return 0.0
-    end
+    q5_2_pow4 = (5 / 2 - q)^4
+    q3_2_pow4 = (3 / 2 - q)^4
+    q1_2_pow4 = (1 / 2 - q)^4
 
-    result = (5 / 2 - q)^4
+    # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl.
+    result = q5_2_pow4
+    result = result - 5 * (q < 3 / 2) * q3_2_pow4
+    result = result + 10 * (q < 1 / 2) * q1_2_pow4
 
-    if q < 3 / 2
-        result -= 5 * (3 / 2 - q)^4
+    # Zero out result if q >= 5/2
+    result = ifelse(q < 5 / 2, normalization_factor(kernel, h) * result, zero(result))
 
-        if q < 1 / 2
-            result += 10 * (1 / 2 - q)^4
-        end
-    end
-
-    return normalization_factor(kernel, h) * result
+    return result
 end
 
 function kernel_deriv(kernel::SchoenbergQuarticSplineKernel, r::Real, h)