@@ -16,138 +16,194 @@ The following smoothing kernels are currently available:
1616
1717Any Kernel with a stability rating of more than '+++' doesn't suffer from pairing-instability.
1818
19- We recommend to use the [ ` WendlandC2Kernel ` ] ( @ref ) for most applications.
19+ We recommend using the [ ` WendlandC2Kernel ` ] ( @ref ) for most applications.
2020If less smoothing is needed, try [ ` SchoenbergCubicSplineKernel ` ] ( @ref ) , for more smoothing try [ ` WendlandC6Kernel ` ] ( @ref ) .
2121
2222``` @eval
23-
24- using TrixiParticles
25- using CairoMakie
26-
27- # --- Group the kernels for combined plotting ---
28- wendland_kernels = [
29- ("Wendland C2", WendlandC2Kernel{2}()),
30- ("Wendland C4", WendlandC4Kernel{2}()),
31- ("Wendland C6", WendlandC6Kernel{2}()),
32- ]
33-
34- schoenberg_kernels = [
35- ("Cubic Spline", SchoenbergCubicSplineKernel{2}()),
36- ("Quartic Spline", SchoenbergQuarticSplineKernel{2}()),
37- ("Quintic Spline", SchoenbergQuinticSplineKernel{2}()),
38- ]
39-
40- other_kernels = [
41- ("Gaussian", GaussianKernel{2}()),
42- ("Poly6", Poly6Kernel{2}()),
43- ("Laguerre-Gauss", LaguerreGaussKernel{2}()),
44- ]
45-
46- spiky_kernel_group = [
47- ("Spiky Kernel", SpikyKernel{2}()),
48- ]
49-
50- # A list of all kernel groups to be plotted
51- # A boolean flag controls whether to apply the consistent y-range
52- kernel_groups = [
53- (title="Wendland Kernels", kernels=wendland_kernels, use_consistent_range=true),
54- (title="Schoenberg Spline Kernels", kernels=schoenberg_kernels, use_consistent_range=true),
55- (title="Other Kernels", kernels=other_kernels, use_consistent_range=true),
56- (title="Spiky Kernel", kernels=spiky_kernel_group, use_consistent_range=false),
57- ]
58-
59- # --- Pre-calculate global y-ranges for consistency ---
60- kernels_for_range_calc = vcat(wendland_kernels, schoenberg_kernels, other_kernels)
61-
62- q_range = range(0, 3, length=300)
63- h = 1.0
64- min_val, max_val = Inf, -Inf
65- min_deriv, max_deriv = Inf, -Inf
66-
67- for (_, kernel_obj) in kernels_for_range_calc
68- kernel_values = [TrixiParticles.kernel(kernel_obj, q, h) for q in q_range]
23+ using TrixiParticles
24+ using CairoMakie
25+
26+ # --- Group the kernels for combined plotting ---
27+ wendland_kernels = [
28+ ("Wendland C2", WendlandC2Kernel{2}()),
29+ ("Wendland C4", WendlandC4Kernel{2}()),
30+ ("Wendland C6", WendlandC6Kernel{2}()),
31+ ]
32+
33+ schoenberg_kernels = [
34+ ("Cubic Spline", SchoenbergCubicSplineKernel{2}()),
35+ ("Quartic Spline", SchoenbergQuarticSplineKernel{2}()),
36+ ("Quintic Spline", SchoenbergQuinticSplineKernel{2}()),
37+ ]
38+
39+ other_kernels = [
40+ ("Gaussian", GaussianKernel{2}()),
41+ ("Poly6", Poly6Kernel{2}()),
42+ ("Laguerre-Gauss", LaguerreGaussKernel{2}()),
43+ ]
44+
45+ spiky_kernel_group = [
46+ ("Spiky Kernel", SpikyKernel{2}()),
47+ ]
48+
49+ # A list of all kernel groups to be plotted
50+ # A boolean flag controls whether to apply the consistent y-range
51+ kernel_groups = [
52+ (title="Wendland Kernels", kernels=wendland_kernels, use_consistent_range=true),
53+ (title="Schoenberg Spline Kernels", kernels=schoenberg_kernels, use_consistent_range=true),
54+ (title="Other Kernels", kernels=other_kernels, use_consistent_range=true),
55+ (title="Spiky Kernel", kernels=spiky_kernel_group, use_consistent_range=false),
56+ ]
57+
58+ # --- Pre-calculate global y-ranges for consistency ---
59+ kernels_for_range_calc = vcat(wendland_kernels, schoenberg_kernels, other_kernels)
60+
61+ q_range = range(0, 3, length=300)
62+ h = 1.0
63+ min_val, max_val = Inf, -Inf
64+ min_deriv, max_deriv = Inf, -Inf
65+
66+ for (_, kernel_obj) in kernels_for_range_calc
67+ kernel_values = [TrixiParticles.kernel(kernel_obj, q, h) for q in q_range]
6968 kernel_derivs = [TrixiParticles.kernel_deriv(kernel_obj, q, h) for q in q_range]
7069
71- global min_val = min(min_val, minimum(kernel_values))
72- global max_val = max(max_val, maximum(kernel_values))
73- global min_deriv = min(min_deriv, minimum(kernel_derivs))
74- global max_deriv = max(max_deriv, maximum(kernel_derivs))
75- end
70+ global min_val = min(min_val, minimum(kernel_values))
71+ global max_val = max(max_val, maximum(kernel_values))
72+ global min_deriv = min(min_deriv, minimum(kernel_derivs))
73+ global max_deriv = max(max_deriv, maximum(kernel_derivs))
74+ end
7675
77- # Add 10% padding to the y-limits for better visuals
78- y_range_val = (min_val - 0.1 * (max_val - min_val),
79- max_val + 0.1 * (max_val - min_val))
80- y_range_deriv = (min_deriv - 0.1 * (max_deriv - min_deriv),
81- max_deriv + 0.1 * (max_deriv - min_deriv))
76+ # Add 10% padding to the y-limits for better visuals
77+ y_range_val = (min_val - 0.1 * (max_val - min_val),
78+ max_val + 0.1 * (max_val - min_val))
79+ y_range_deriv = (min_deriv - 0.1 * (max_deriv - min_deriv),
80+ max_deriv + 0.1 * (max_deriv - min_deriv))
8281
83- fig = Figure(size = (1000, 1200), fontsize=16)
82+ fig = Figure(size = (1000, 1200), fontsize=16)
8483
85- for (i, group) in enumerate(kernel_groups)
86- ax_val = Axis(fig[i, 1],
87- xlabel = "q = r/h", ylabel = "w(q)",
88- title = group.title)
84+ for (i, group) in enumerate(kernel_groups)
85+ ax_val = Axis(fig[i, 1],
86+ xlabel = "q = r/h", ylabel = "w(q)",
87+ title = group.title)
8988
90- ax_deriv = Axis(fig[i, 2],
91- xlabel = "q = r/h", ylabel = "w'(q)",
92- title = group.title)
89+ ax_deriv = Axis(fig[i, 2],
90+ xlabel = "q = r/h", ylabel = "w'(q)",
91+ title = group.title)
9392
94- if group.use_consistent_range
95- ylims!(ax_val, y_range_val)
96- ylims!(ax_deriv, y_range_deriv)
97- end
93+ if group.use_consistent_range
94+ ylims!(ax_val, y_range_val)
95+ ylims!(ax_deriv, y_range_deriv)
96+ end
9897
99- hlines!(ax_val, [0.0], linestyle = :dash)
100- hlines!(ax_deriv, [0.0], linestyle = :dash)
98+ hlines!(ax_val, [0.0], linestyle = :dash)
99+ hlines!(ax_deriv, [0.0], linestyle = :dash)
101100
102- for (name, kernel_obj) in group.kernels
103- kernel_values = [TrixiParticles.kernel(kernel_obj, q, h) for q in q_range]
104- kernel_derivs = [TrixiParticles.kernel_deriv(kernel_obj, q, h) for q in q_range]
101+ for (name, kernel_obj) in group.kernels
102+ kernel_values = [TrixiParticles.kernel(kernel_obj, q, h) for q in q_range]
103+ kernel_derivs = [TrixiParticles.kernel_deriv(kernel_obj, q, h) for q in q_range]
105104
106- lines!(ax_val, q_range, kernel_values, label=name, linewidth=2.5)
107- lines!(ax_deriv, q_range, kernel_derivs, label=name, linewidth=2.5)
108- end
105+ lines!(ax_val, q_range, kernel_values, label=name, linewidth=2.5)
106+ lines!(ax_deriv, q_range, kernel_derivs, label=name, linewidth=2.5)
107+ end
109108
110- axislegend(ax_val, position = :rt)
111- axislegend(ax_deriv, position = :rt)
112- end
109+ axislegend(ax_val, position = :rt)
110+ axislegend(ax_deriv, position = :rt)
111+ end
113112
114- # Add row gaps between the 4 rows (3 gaps total)
115- for i in 1:(length(kernel_groups) - 1)
116- rowgap!(fig.layout, i, 25)
117- end
113+ # Add row gaps between the 4 rows (3 gaps total)
114+ for i in 1:(length(kernel_groups) - 1)
115+ rowgap!(fig.layout, i, 25)
116+ end
118117
119118CairoMakie.save("smoothing_kernels.png", fig)
120-
119+
121120```
122121
123122![ Radial profiles and derivatives of the available smoothing kernels] ( smoothing_kernels.png )
124123
125124
126125!!! note "Usage"
127126 The kernel can be called as
128- ```
127+ ```jldoctest kernels; output = false, setup = :(using TrixiParticles, LinearAlgebra; smoothing_kernel = WendlandC2Kernel{2}(); h = 1.0; r = 0.5)
129128 TrixiParticles.kernel(smoothing_kernel, r, h)
129+
130+ # output
131+ 0.35250333098869013
130132 ```
131133 The length of the compact support can be obtained as
132- ```
134+ ```jldoctest kernels; output = false
133135 TrixiParticles.compact_support(smoothing_kernel, h)
136+
137+ # output
138+ 2.0
134139 ```
135140
136141 Note that ``r`` has to be a scalar, so in the context of SPH, the kernel
137142 should be used as
138143 ```math
139- W(\Vert r_a - r_b \Vert, h).
144+ W(\Vert r_a - r_b \Vert, h).
140145 ```
141146
142147 The gradient required in SPH,
143148 ```math
144149 \nabla_{r_a} W(\Vert r_a - r_b \Vert, h)
145150 ```
146151 can be called as
147- ```
152+ ```jldoctest kernels; output = false, setup = :(pos_diff = SVector(0.5, 0.5); distance = norm(pos_diff))
148153 TrixiParticles.kernel_grad(smoothing_kernel, pos_diff, distance, h)
154+
155+ # output
156+ 2-element SVector{2, Float64} with indices SOneTo(2):
157+ -0.3762063922043116
158+ -0.3762063922043116
159+ ```
160+ where `pos_diff` is $r_a - r_b$ and `distance` is $\Vert r_a - r_b \Vert$,
161+ although in most cases, it is recommended to use
162+ ```jldoctest kernels; output = false, setup=:(trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), sol=nothing); system = fluid_system; particle = 1)
163+ TrixiParticles.smoothing_kernel_grad(system, pos_diff, distance, particle)
164+
165+ # output
166+ 2-element SVector{2, Float64} with indices SOneTo(2):
167+ 0.0
168+ 0.0
169+ ```
170+ instead, which will also take into account if the system is using a corrected kernel.
171+
172+ The equivalent function to obtain the kernel value is
173+ ```jldoctest kernels; output = false
174+ TrixiParticles.smoothing_kernel(system, distance, particle)
175+
176+ # output
177+ 0.0
178+ ```
179+
180+ In performance-critical code, it is recommended to use the `_unsafe` versions of these
181+ functions, which skip the compact support and zero distance checks:
182+ ```jldoctest kernels; output = false
183+ TrixiParticles.smoothing_kernel_unsafe(system, distance, particle)
184+ TrixiParticles.smoothing_kernel_grad_unsafe(system, pos_diff, distance, particle)
185+
186+ # output
187+ 2-element SVector{2, Float64} with indices SOneTo(2):
188+ -106899.65266169122
189+ -106899.65266169122
190+ ```
191+ These functions are only safe to use if compact support and zero distance have already
192+ been checked.
193+ The compact support is automatically checked when inside a `foreach_point_neighbor`
194+ or `foreach_neighbor` loop.
195+ For the zero distance check, use
196+ ```jldoctest kernels; output = false, setup = :(almostzero = 1e-12)
197+ TrixiParticles.skip_zero_distance(system) && distance < almostzero
198+
199+ # output
200+ false
149201 ```
150- where `pos_diff` is $r_a - r_b$ and `distance` is $\Vert r_a - r_b \Vert$.
202+ where `skip_zero_distance` returns `true` if the kernel gradient is zero at zero
203+ distance, which is not the case for some corrected kernels.
204+ Here, `almostzero` is a small threshold, which should be chosen as `sqrt(eps(h^2))`,
205+ where `h` is the smoothing length, since `distance^2` is in the order of `h^2`.
206+ Note that `sqrt(eps(h^2)) != eps(h)`.
151207
152208``` @autodocs
153209Modules = [TrixiParticles]
0 commit comments