Note

Go to the end to download the full example code.

SPH kernels#

This example shows the all the SPH kernels

 8 import matplotlib.pyplot as plt
 9 import numpy as np
10
11 import shamrock

utilities & check integral == 1

17 def compute_integ_3d(q, W):
18     if hasattr(np, "trapezoid"):
19         integrate_func = getattr(np, "trapezoid")
20     else:
21         integrate_func = getattr(np, "trapz")
22     return integrate_func(4 * np.pi * q**2 * W, q)
23
24
25 def plot_test_sph_kernel(q, f, df, W, dW, title, ax):
26
27     integral_result = compute_integ_3d(q, W)
28     assert np.abs(integral_result - 1) < 1e-6, (
29         "3D integration of 4π q² W(q) is not 1, kernel: " + title
30     )
31
32     ax[0].plot(q, W, label=f"$W_{{{title}}}(q)$")
33     ax[1].plot(q, dW, label=f"$dW_{{{title}}}(q)$")
34
35
36 q = np.linspace(0, 4, 1000)

Cubic splines

41 fig, axs = plt.subplots(1, 2, figsize=(10, 5))
42
43 f_M4 = [shamrock.math.sphkernel.M4_f(x) for x in q]
44 w3d_M4 = [shamrock.math.sphkernel.M4_W3d(x, 1) for x in q]
45 df_M4 = [shamrock.math.sphkernel.M4_df(x) for x in q]
46 dW3d_M4 = [shamrock.math.sphkernel.M4_dW3d(x, 1) for x in q]
47 plot_test_sph_kernel(q, f_M4, df_M4, w3d_M4, dW3d_M4, "M4", axs)
48
49 f_M6 = [shamrock.math.sphkernel.M6_f(x) for x in q]
50 w3d_M6 = [shamrock.math.sphkernel.M6_W3d(x, 1) for x in q]
51 df_M6 = [shamrock.math.sphkernel.M6_df(x) for x in q]
52 dW3d_M6 = [shamrock.math.sphkernel.M6_dW3d(x, 1) for x in q]
53 plot_test_sph_kernel(q, f_M6, df_M6, w3d_M6, dW3d_M6, "M6", axs)
54
55 f_M8 = [shamrock.math.sphkernel.M8_f(x) for x in q]
56 w3d_M8 = [shamrock.math.sphkernel.M8_W3d(x, 1) for x in q]
57 df_M8 = [shamrock.math.sphkernel.M8_df(x) for x in q]
58 dW3d_M8 = [shamrock.math.sphkernel.M8_dW3d(x, 1) for x in q]
59 plot_test_sph_kernel(q, f_M8, df_M8, w3d_M8, dW3d_M8, "M8", axs)
60
61 axs[0].legend()
62 axs[1].legend()
63 axs[0].set_xlabel(r"$q$")
64 axs[1].set_xlabel(r"$q$")
65 plt.show()
run show all sph kernels

Wendland kernels

71 fig, axs = plt.subplots(1, 2, figsize=(10, 5))
72
73 f_C2 = [shamrock.math.sphkernel.C2_f(x) for x in q]
74 w3d_C2 = [shamrock.math.sphkernel.C2_W3d(x, 1) for x in q]
75 df_C2 = [shamrock.math.sphkernel.C2_df(x) for x in q]
76 dW3d_C2 = [shamrock.math.sphkernel.C2_dW3d(x, 1) for x in q]
77 plot_test_sph_kernel(q, f_C2, df_C2, w3d_C2, dW3d_C2, "C2", axs)
78
79 f_C4 = [shamrock.math.sphkernel.C4_f(x) for x in q]
80 w3d_C4 = [shamrock.math.sphkernel.C4_W3d(x, 1) for x in q]
81 df_C4 = [shamrock.math.sphkernel.C4_df(x) for x in q]
82 dW3d_C4 = [shamrock.math.sphkernel.C4_dW3d(x, 1) for x in q]
83 plot_test_sph_kernel(q, f_C4, df_C4, w3d_C4, dW3d_C4, "C4", axs)
84
85 f_C6 = [shamrock.math.sphkernel.C6_f(x) for x in q]
86 w3d_C6 = [shamrock.math.sphkernel.C6_W3d(x, 1) for x in q]
87 df_C6 = [shamrock.math.sphkernel.C6_df(x) for x in q]
88 dW3d_C6 = [shamrock.math.sphkernel.C6_dW3d(x, 1) for x in q]
89 plot_test_sph_kernel(q, f_C6, df_C6, w3d_C6, dW3d_C6, "C6", axs)
90
91 axs[0].legend()
92 axs[1].legend()
93 axs[0].set_xlabel(r"$q$")
94 axs[1].set_xlabel(r"$q$")
95 plt.show()
run show all sph kernels

Double hump kernels

101 fig, axs = plt.subplots(1, 2, figsize=(10, 5))
102
103 f_M4DH = [shamrock.math.sphkernel.M4DH_f(x) for x in q]
104 w3d_M4DH = [shamrock.math.sphkernel.M4DH_W3d(x, 1) for x in q]
105 df_M4DH = [shamrock.math.sphkernel.M4DH_df(x) for x in q]
106 dW3d_M4DH = [shamrock.math.sphkernel.M4DH_dW3d(x, 1) for x in q]
107 plot_test_sph_kernel(q, f_M4DH, df_M4DH, w3d_M4DH, dW3d_M4DH, "M4DH", axs)
108
109 f_M4DH3 = [shamrock.math.sphkernel.M4DH3_f(x) for x in q]
110 w3d_M4DH3 = [shamrock.math.sphkernel.M4DH3_W3d(x, 1) for x in q]
111 df_M4DH3 = [shamrock.math.sphkernel.M4DH3_df(x) for x in q]
112 dW3d_M4DH3 = [shamrock.math.sphkernel.M4DH3_dW3d(x, 1) for x in q]
113 plot_test_sph_kernel(q, f_M4DH3, df_M4DH3, w3d_M4DH3, dW3d_M4DH3, "M4DH3", axs)
114
115 f_M4DH5 = [shamrock.math.sphkernel.M4DH5_f(x) for x in q]
116 w3d_M4DH5 = [shamrock.math.sphkernel.M4DH5_W3d(x, 1) for x in q]
117 df_M4DH5 = [shamrock.math.sphkernel.M4DH5_df(x) for x in q]
118 dW3d_M4DH5 = [shamrock.math.sphkernel.M4DH5_dW3d(x, 1) for x in q]
119 plot_test_sph_kernel(q, f_M4DH5, df_M4DH5, w3d_M4DH5, dW3d_M4DH5, "M4DH5", axs)
120
121 f_M4DH7 = [shamrock.math.sphkernel.M4DH7_f(x) for x in q]
122 w3d_M4DH7 = [shamrock.math.sphkernel.M4DH7_W3d(x, 1) for x in q]
123 df_M4DH7 = [shamrock.math.sphkernel.M4DH7_df(x) for x in q]
124 dW3d_M4DH7 = [shamrock.math.sphkernel.M4DH7_dW3d(x, 1) for x in q]
125 plot_test_sph_kernel(q, f_M4DH7, df_M4DH7, w3d_M4DH7, dW3d_M4DH7, "M4DH7", axs)
126
127 axs[0].legend()
128 axs[1].legend()
129 axs[0].set_xlabel(r"$q$")
130 axs[1].set_xlabel(r"$q$")
131 plt.show()
run show all sph kernels

shifted cubic splines kernels

137 fig, axs = plt.subplots(1, 2, figsize=(10, 5))
138
139 f_M4Shift2 = [shamrock.math.sphkernel.M4Shift2_f(x) for x in q]
140 w3d_M4Shift2 = [shamrock.math.sphkernel.M4Shift2_W3d(x, 1) for x in q]
141 df_M4Shift2 = [shamrock.math.sphkernel.M4Shift2_df(x) for x in q]
142 dW3d_M4Shift2 = [shamrock.math.sphkernel.M4Shift2_dW3d(x, 1) for x in q]
143 plot_test_sph_kernel(q, f_M4Shift2, df_M4Shift2, w3d_M4Shift2, dW3d_M4Shift2, "M4Shift2", axs)
144
145 f_M4Shift4 = [shamrock.math.sphkernel.M4Shift4_f(x) for x in q]
146 w3d_M4Shift4 = [shamrock.math.sphkernel.M4Shift4_W3d(x, 1) for x in q]
147 df_M4Shift4 = [shamrock.math.sphkernel.M4Shift4_df(x) for x in q]
148 dW3d_M4Shift4 = [shamrock.math.sphkernel.M4Shift4_dW3d(x, 1) for x in q]
149 plot_test_sph_kernel(q, f_M4Shift4, df_M4Shift4, w3d_M4Shift4, dW3d_M4Shift4, "M4Shift4", axs)
150
151 f_M4Shift8 = [shamrock.math.sphkernel.M4Shift8_f(x) for x in q]
152 w3d_M4Shift8 = [shamrock.math.sphkernel.M4Shift8_W3d(x, 1) for x in q]
153 df_M4Shift8 = [shamrock.math.sphkernel.M4Shift8_df(x) for x in q]
154 dW3d_M4Shift8 = [shamrock.math.sphkernel.M4Shift8_dW3d(x, 1) for x in q]
155 plot_test_sph_kernel(q, f_M4Shift8, df_M4Shift8, w3d_M4Shift8, dW3d_M4Shift8, "M4Shift8", axs)
156
157 f_M4Shift16 = [shamrock.math.sphkernel.M4Shift16_f(x) for x in q]
158 w3d_M4Shift16 = [shamrock.math.sphkernel.M4Shift16_W3d(x, 1) for x in q]
159 df_M4Shift16 = [shamrock.math.sphkernel.M4Shift16_df(x) for x in q]
160 dW3d_M4Shift16 = [shamrock.math.sphkernel.M4Shift16_dW3d(x, 1) for x in q]
161 plot_test_sph_kernel(q, f_M4Shift16, df_M4Shift16, w3d_M4Shift16, dW3d_M4Shift16, "M4Shift16", axs)
162
163 axs[0].legend()
164 axs[1].legend()
165 axs[0].set_xlabel(r"$q$")
166 axs[1].set_xlabel(r"$q$")
167 plt.show()
run show all sph kernels

Total running time of the script: (0 minutes 2.192 seconds)

Estimated memory usage: 108 MB

Gallery generated by Sphinx-Gallery