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

Use shamrock documentation style for matplotlib

15 shamrock.matplotlib.set_shamrock_mpl_style()

utilities & check integral == 1

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

Cubic splines

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

Wendland kernels

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

Double hump kernels

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

shifted cubic splines kernels

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

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

Estimated memory usage: 159 MB

Gallery generated by Sphinx-Gallery