Note
Go to the end to download the full example code.
FMM math demo#
This example shows how to use the FMM maths to compute the force between two points
As always, we start by importing the necessary libraries
10 import matplotlib
11 import matplotlib.pyplot as plt
12 import numpy as np
13 from mpl_toolkits.mplot3d.art3d import Line3DCollection, Poly3DCollection
14
15 import shamrock
Utilities#
You can ignore this first block, it just contains some utility functions to draw the AABB and the arrows We only defines the function draw_aabb and draw_arrow, which are used to draw the AABB and the arrows in the plots and the function draw_box_pair, which is used to draw the box pair with all the vectors needed to compute the FMM force
Click here to expand the utility code
32 def draw_aabb(ax, aabb, color, alpha):
33 """
34 Draw a 3D AABB in matplotlib
35
36 Parameters
37 ----------
38 ax : matplotlib.Axes3D
39 The axis to draw the AABB on
40 aabb : shamrock.math.AABB_f64_3
41 The AABB to draw
42 color : str
43 The color of the AABB
44 alpha : float
45 The transparency of the AABB
46 """
47 xmin, ymin, zmin = aabb.lower
48 xmax, ymax, zmax = aabb.upper
49
50 points = [
51 aabb.lower,
52 (aabb.lower[0], aabb.lower[1], aabb.upper[2]),
53 (aabb.lower[0], aabb.upper[1], aabb.lower[2]),
54 (aabb.lower[0], aabb.upper[1], aabb.upper[2]),
55 (aabb.upper[0], aabb.lower[1], aabb.lower[2]),
56 (aabb.upper[0], aabb.lower[1], aabb.upper[2]),
57 (aabb.upper[0], aabb.upper[1], aabb.lower[2]),
58 aabb.upper,
59 ]
60
61 faces = [
62 [points[0], points[1], points[3], points[2]],
63 [points[4], points[5], points[7], points[6]],
64 [points[0], points[1], points[5], points[4]],
65 [points[2], points[3], points[7], points[6]],
66 [points[0], points[2], points[6], points[4]],
67 [points[1], points[3], points[7], points[5]],
68 ]
69
70 edges = [
71 [points[0], points[1]],
72 [points[0], points[2]],
73 [points[0], points[4]],
74 [points[1], points[3]],
75 [points[1], points[5]],
76 [points[2], points[3]],
77 [points[2], points[6]],
78 [points[3], points[7]],
79 [points[4], points[5]],
80 [points[4], points[6]],
81 [points[5], points[7]],
82 [points[6], points[7]],
83 ]
84
85 collection = Poly3DCollection(faces, alpha=alpha, color=color)
86 ax.add_collection3d(collection)
87
88 edge_collection = Line3DCollection(edges, color="k", alpha=alpha)
89 ax.add_collection3d(edge_collection)
90
91
92 def draw_arrow(ax, p1, p2, color, label, arr_scale=0.1):
93 length = np.linalg.norm(np.array(p2) - np.array(p1))
94 arrow_length_ratio = arr_scale / length
95 ax.quiver(
96 p1[0],
97 p1[1],
98 p1[2],
99 p2[0] - p1[0],
100 p2[1] - p1[1],
101 p2[2] - p1[2],
102 color=color,
103 label=label,
104 arrow_length_ratio=arrow_length_ratio,
105 )
FMM force computation#
Let’s start by assuming that we have two particles at positions \(\mathbf{x}_i\) and \(\mathbf{x}_j\) contained in two boxes (\(A\) and \(B\)) whose centers are at positions \(\mathbf{s}_a\) and \(\mathbf{s}_b\) respectively. The positions of the particles relative to their respective boxes are then:
\[\begin{split}\mathbf{a}_i = \mathbf{x}_i - \mathbf{s}_a \\
\mathbf{b}_j = \mathbf{x}_j - \mathbf{s}_b\end{split}\]
and the distance between the centers of the boxes is:
\[\mathbf{r} = \mathbf{s}_b - \mathbf{s}_a\]
This implies that the distance between the two particles is:
\[\mathbf{x}_j - \mathbf{x}_i = \mathbf{r} + \mathbf{b}_j - \mathbf{a}_i\]
If we denote the Green function for an inverse distance \(G(\mathbf{x}) = 1 / \vert\vert\mathbf{x}\vert\vert\), then the potential exerted onto particle \(i\) is:
\[\begin{split}\Phi_i = \Phi (\mathbf{x}_i) &= \int - \frac{\mathcal{G} \rho(\mathbf{x}_j)}{\vert\vert\mathbf{x}_i - \mathbf{x}_j\vert\vert} d\mathbf{x}_j \\
&= - \mathcal{G} \int \rho(\mathbf{x}_j) G(\mathbf{x}_j - \mathbf{x}_i) d\mathbf{x}_j\end{split}\]
and the force exerted onto particle \(i\) is:
\[\begin{split}\mathbf{f}_i = -\nabla_i \Phi (\mathbf{x}_i) &= \int - \nabla_i \frac{\mathcal{G} \rho(\mathbf{x}_j)}{\vert\vert\mathbf{x}_i - \mathbf{x}_j\vert\vert} d\mathbf{x}_j \\
&= \mathcal{G} \int \rho(\mathbf{x}_j) \nabla_i G(\mathbf{x}_j - \mathbf{x}_i) d\mathbf{x}_j \\
&= -\mathcal{G} \int \rho(\mathbf{x}_j) \nabla_j G(\mathbf{x}_j - \mathbf{x}_i) d\mathbf{x}_j\end{split}\]
Now let’s expand the green function in a Taylor series to order \(p\).
\[\begin{split}G(\mathbf{x}_j - \mathbf{x}_i) &= G(\mathbf{r} + \mathbf{b}_j - \mathbf{a}_i) \\
&= \sum_{k = 0}^p \frac{(-1)^k}{k!} \mathbf{a}_i^{(k)} \cdot \sum_{n=0}^{p-k} \frac{1}{n!} D_{n+k} \cdot \mathbf{b}_j^{(n)}\end{split}\]
where \(D_{n} = \nabla^{(n)}_r G(\mathbf{r})\) is the n-th order derivative of the Green function and the operator \(\mathbf{a}_i^{(k)}\) is the tensor product of \(k\) \(\mathbf{a}_i\) vectors.
Similarly the gradient of the green function is:
\[\begin{split}\nabla_j G(\mathbf{x}_j - \mathbf{x}_i) &= \nabla_r G(\mathbf{r} + \mathbf{b}_j - \mathbf{a}_i) \\
&= \sum_{k = 0}^p \frac{(-1)^k}{k!} \mathbf{a}_i^{(k)} \cdot \sum_{n=0}^{p-k} \frac{1}{n!} D_{n+k+1} \cdot \mathbf{b}_j^{(n)}\end{split}\]
Now we can plug that back into the expression for the force & potential:
\[\begin{split}\Phi_i &= - \mathcal{G} \int \rho(\mathbf{x}_j) G(\mathbf{x}_j - \mathbf{x}_i) d\mathbf{x}_j \\
&= - \mathcal{G} \sum_{k = 0}^p \frac{1}{k!} \mathbf{a}_i^{(k)} \cdot \underbrace{(-1)^k \sum_{n=0}^{p-k} \frac{1}{n!} D_{n+k} \cdot \underbrace{\int \rho(\mathbf{x}_j) \mathbf{b}_j^{(n)} d\mathbf{x}_j}_{Q_n^B}}_{M_{p,k}} \\\end{split}\]
\[\begin{split}\mathbf{f}_i &= -\mathcal{G} \int \rho(\mathbf{x}_j) \nabla_j G(\mathbf{x}_j - \mathbf{x}_i) d\mathbf{x}_j \\
&= -\mathcal{G} \sum_{k = 0}^p \frac{1}{k!} \mathbf{a}_i^{(k)} \cdot \underbrace{(-1)^k \sum_{n=0}^{p-k} \frac{1}{n!} D_{n+k+1} \cdot \underbrace{\int \rho(\mathbf{x}_j)\mathbf{b}_j^{(n)} d\mathbf{x}_j}_{Q_n^B}}_{dM_{p,k} = M_{p+1,k+1}}\end{split}\]
As one can tell sadly the two expressions while similar do not share the same terms.
I will not go in this rabit hole of using the same expansion for both now but the idea is to use the primitive of the force which is the same expansion as the force but with the primitive of \(\mathbf{a}_i^{(k)}\) instead.
\[\Phi_i = - \int \mathbf{f}_i = -\mathcal{G} \sum_{k = 0}^p \frac{1}{k!} \int\mathbf{a}_i^{(k)} \cdot {M_{p+1,k+1}}\]
Mass moments#
def plot_mass_moment_case(s_B,box_B_size,x_j):
204 def plot_mass_moment_case(s_B, box_B_size, x_j):
205 box_B = shamrock.math.AABB_f64_3(
206 (
207 s_B[0] - box_B_size / 2,
208 s_B[1] - box_B_size / 2,
209 s_B[2] - box_B_size / 2,
210 ),
211 (
212 s_B[0] + box_B_size / 2,
213 s_B[1] + box_B_size / 2,
214 s_B[2] + box_B_size / 2,
215 ),
216 )
217
218 fig = plt.figure()
219 ax = fig.add_subplot(111, projection="3d")
220
221 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
222
223 ax.scatter(s_B[0], s_B[1], s_B[2], color="black", label="s_B")
224
225 ax.scatter(x_j[0], x_j[1], x_j[2], color="red", label="$x_j$")
226
227 draw_aabb(ax, box_B, "blue", 0.2)
228
229 center_view = (0.0, 0.0, 0.0)
230 view_size = 2.0
231 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
232 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
233 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
234 ax.set_xlabel("X")
235 ax.set_ylabel("Y")
236 ax.set_zlabel("Z")
237 return ax
Let’s start with the following
248 s_B = (0, 0, 0)
249 box_B_size = 1
250
251 x_j = (0.2, 0.2, 0.2)
252 m_j = 1
253
254 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
255
256 ax = plot_mass_moment_case(s_B, box_B_size, x_j)
257 plt.title("Mass moment illustration")
258 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
259 plt.show()
Here the mass moment of a set of particles (here only one) of mass \(m_j\) is
\[\begin{split}{Q_n^B} &= \int \rho(\mathbf{x}_j) \mathbf{b}_j^{(n)} d\mathbf{x}_j\\
&= \sum_j m_j \mathbf{b}_j^{(n)}\end{split}\]
In Shamrock python bindings the function
will return the collection of symetrical tensors \(\mathbf{b}_j^{(n)}\) for n in between <low order> and <high order> Here are the values of the tensors \({Q_n^B}\) from order 0 up to 5 using shamrock symmetrical tensor collections
Q_n_B = SymTensorCollection_f64_0_5(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0.2),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0.04000000000000001, v_11=0.04000000000000001, v_12=0.04000000000000001, v_22=0.04000000000000001),
t3=SymTensor3d_3(v_000=0.008000000000000002, v_001=0.008000000000000002, v_002=0.008000000000000002, v_011=0.008000000000000002, v_012=0.008000000000000002, v_022=0.008000000000000002, v_111=0.008000000000000002, v_112=0.008000000000000002, v_122=0.008000000000000002, v_222=0.008000000000000002),
t4=SymTensor3d_4(v_0000=0.0016000000000000005, v_0001=0.0016000000000000005, v_0002=0.0016000000000000005, v_0011=0.0016000000000000005, v_0012=0.0016000000000000005, v_0022=0.0016000000000000005, v_0111=0.0016000000000000005, v_0112=0.0016000000000000005, v_0122=0.0016000000000000005, v_0222=0.0016000000000000005, v_1111=0.0016000000000000005, v_1112=0.0016000000000000005, v_1122=0.0016000000000000005, v_1222=0.0016000000000000005, v_2222=0.0016000000000000005),
t5=SymTensor3d_5(v_00000=0.00032000000000000013, v_00001=0.00032000000000000013, v_00002=0.00032000000000000013, v_00011=0.00032000000000000013, v_00012=0.00032000000000000013, v_00022=0.00032000000000000013, v_00111=0.00032000000000000013, v_00112=0.00032000000000000013, v_00122=0.00032000000000000013, v_00222=0.00032000000000000013, v_01111=0.00032000000000000013, v_01112=0.00032000000000000013, v_01122=0.00032000000000000013, v_01222=0.00032000000000000013, v_02222=0.00032000000000000013, v_11111=0.00032000000000000013, v_11112=0.00032000000000000013, v_11122=0.00032000000000000013, v_11222=0.00032000000000000013, v_12222=0.00032000000000000013, v_22222=0.00032000000000000013)
)
Now if we take a displacment that is only along the x axis we get null components in the Q_n_B if for cases that do not only exhibit x
288 s_B = (0, 0, 0)
289 box_B_size = 1
290
291 x_j = (0.5, 0.0, 0.0)
292 m_j = 1
293
294 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
295
296 ax = plot_mass_moment_case(s_B, box_B_size, x_j)
297 plt.title("Mass moment illustration")
298 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
299 plt.show()
300
301 Q_n_B = shamrock.math.SymTensorCollection_f64_0_5.from_vec(b_j)
302 Q_n_B *= m_j
303
304 print("Q_n_B =", Q_n_B)
Q_n_B = SymTensorCollection_f64_0_5(
t0=1,
t1=SymTensor3d_1(v_0=0.5, v_1=0, v_2=0),
t2=SymTensor3d_2(v_00=0.25, v_01=0, v_02=0, v_11=0, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0.125, v_001=0, v_002=0, v_011=0, v_012=0, v_022=0, v_111=0, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=0.0625, v_0001=0, v_0002=0, v_0011=0, v_0012=0, v_0022=0, v_0111=0, v_0112=0, v_0122=0, v_0222=0, v_1111=0, v_1112=0, v_1122=0, v_1222=0, v_2222=0),
t5=SymTensor3d_5(v_00000=0.03125, v_00001=0, v_00002=0, v_00011=0, v_00012=0, v_00022=0, v_00111=0, v_00112=0, v_00122=0, v_00222=0, v_01111=0, v_01112=0, v_01122=0, v_01222=0, v_02222=0, v_11111=0, v_11112=0, v_11122=0, v_11222=0, v_12222=0, v_22222=0)
)
Gravitational moments#
def plot_mass_moment_case(s_B,box_B_size,x_j):
319 def plot_grav_moment_case(s_A, box_A_size, s_B, box_B_size, x_j):
320 box_A = shamrock.math.AABB_f64_3(
321 (
322 s_A[0] - box_A_size / 2,
323 s_A[1] - box_A_size / 2,
324 s_A[2] - box_A_size / 2,
325 ),
326 (
327 s_A[0] + box_A_size / 2,
328 s_A[1] + box_A_size / 2,
329 s_A[2] + box_A_size / 2,
330 ),
331 )
332
333 box_B = shamrock.math.AABB_f64_3(
334 (
335 s_B[0] - box_B_size / 2,
336 s_B[1] - box_B_size / 2,
337 s_B[2] - box_B_size / 2,
338 ),
339 (
340 s_B[0] + box_B_size / 2,
341 s_B[1] + box_B_size / 2,
342 s_B[2] + box_B_size / 2,
343 ),
344 )
345
346 fig = plt.figure()
347 ax = fig.add_subplot(111, projection="3d")
348
349 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
350
351 draw_arrow(ax, s_A, s_B, "purple", "$r = s_B - s_A$")
352
353 ax.scatter(s_A[0], s_A[1], s_A[2], color="black", label="s_A")
354
355 ax.scatter(s_B[0], s_B[1], s_B[2], color="green", label="s_B")
356
357 ax.scatter(x_j[0], x_j[1], x_j[2], color="red", label="$x_j$")
358
359 draw_aabb(ax, box_A, "blue", 0.1)
360 draw_aabb(ax, box_B, "red", 0.1)
361
362 center_view = (0.5, 0.0, 0.0)
363 view_size = 2.0
364 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
365 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
366 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
367 ax.set_xlabel("X")
368 ax.set_ylabel("Y")
369 ax.set_zlabel("Z")
370 return ax
Let’s now show the example of a gravitational moment, for the following case
380 s_B = (0, 0, 0)
381 s_A = (1, 0, 0)
382
383 box_B_size = 0.5
384 box_A_size = 0.5
385
386 x_j = (0.2, 0.2, 0.0)
387 m_j = 1
388
389 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
390 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
391
392 ax = plot_grav_moment_case(s_A, box_A_size, s_B, box_B_size, x_j)
393 plt.title("Grav moment illustration")
394 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
395 plt.show()
The mass moment \({Q_n^B}\) is
Q_n_B = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0.008000000000000002, v_001=0.008000000000000002, v_002=0, v_011=0.008000000000000002, v_012=0, v_022=0, v_111=0.008000000000000002, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=0.0016000000000000005, v_0001=0.0016000000000000005, v_0002=0, v_0011=0.0016000000000000005, v_0012=0, v_0022=0, v_0111=0.0016000000000000005, v_0112=0, v_0122=0, v_0222=0, v_1111=0.0016000000000000005, v_1112=0, v_1122=0, v_1222=0, v_2222=0)
)
The green function n’th gradients \(D_{n+k+1}\) are
405 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
406 print("D_n =", D_n)
D_n = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=1, v_1=-0, v_2=-0),
t2=SymTensor3d_2(v_00=2, v_01=-0, v_02=-0, v_11=-1, v_12=0, v_22=-1),
t3=SymTensor3d_3(v_000=6, v_001=0, v_002=0, v_011=-3, v_012=0, v_022=-3, v_111=0, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=24, v_0001=0, v_0002=0, v_0011=-12, v_0012=0, v_0022=-12, v_0111=0, v_0112=0, v_0122=0, v_0222=0, v_1111=9, v_1112=0, v_1122=3, v_1222=0, v_2222=9),
t5=SymTensor3d_5(v_00000=120, v_00001=-0, v_00002=-0, v_00011=-60, v_00012=0, v_00022=-60, v_00111=0, v_00112=0, v_00122=0, v_00222=0, v_01111=45, v_01112=-0, v_01122=15, v_01222=-0, v_02222=45, v_11111=0, v_11112=-0, v_11122=0, v_11222=0, v_12222=-0, v_22222=0)
)
And finally the gravitational moments \(dM_{p,k}\) are
410 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
411 print("dM_k =", dM_k)
dM_k = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=1.431, v_1=-0.3600000000000001, v_2=0),
t2=SymTensor3d_2(v_00=-3.3600000000000003, v_01=1.26, v_02=0, v_11=1.56, v_12=-0, v_22=1.8000000000000003),
t3=SymTensor3d_3(v_000=12, v_001=-4.800000000000001, v_002=0, v_011=-5.7, v_012=0, v_022=-6.300000000000001, v_111=3.6000000000000005, v_112=0, v_122=1.2000000000000002, v_222=0),
t4=SymTensor3d_4(v_0000=-48, v_0001=12, v_0002=-0, v_0011=24, v_0012=-0, v_0022=24, v_0111=-9, v_0112=-0, v_0122=-3, v_0222=-0, v_1111=-18, v_1112=-0, v_1122=-6, v_1222=-0, v_2222=-18),
t5=SymTensor3d_5(v_00000=120, v_00001=-0, v_00002=-0, v_00011=-60, v_00012=0, v_00022=-60, v_00111=0, v_00112=0, v_00122=0, v_00222=0, v_01111=45, v_01112=-0, v_01122=15, v_01222=-0, v_02222=45, v_11111=0, v_11112=-0, v_11122=0, v_11222=0, v_12222=-0, v_22222=0)
)
From Gravitational moments to force#
def plot_fmm_case(s_A,box_A_size,x_i,s_B,box_B_size,x_j, f_i_fmm, f_i_exact):
425 def plot_fmm_case(s_A, box_A_size, x_i, s_B, box_B_size, x_j, f_i_fmm, f_i_exact, fscale_fact):
426 box_A = shamrock.math.AABB_f64_3(
427 (
428 s_A[0] - box_A_size / 2,
429 s_A[1] - box_A_size / 2,
430 s_A[2] - box_A_size / 2,
431 ),
432 (
433 s_A[0] + box_A_size / 2,
434 s_A[1] + box_A_size / 2,
435 s_A[2] + box_A_size / 2,
436 ),
437 )
438
439 box_B = shamrock.math.AABB_f64_3(
440 (
441 s_B[0] - box_B_size / 2,
442 s_B[1] - box_B_size / 2,
443 s_B[2] - box_B_size / 2,
444 ),
445 (
446 s_B[0] + box_B_size / 2,
447 s_B[1] + box_B_size / 2,
448 s_B[2] + box_B_size / 2,
449 ),
450 )
451
452 fig = plt.figure()
453 ax = fig.add_subplot(111, projection="3d")
454
455 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
456 draw_arrow(ax, s_A, x_i, "blue", "$a_i = x_i - s_A$")
457
458 draw_arrow(ax, s_A, s_B, "purple", "$r = s_B - s_A$")
459
460 ax.scatter(s_A[0], s_A[1], s_A[2], color="black", label="s_A")
461
462 ax.scatter(s_B[0], s_B[1], s_B[2], color="green", label="s_B")
463
464 ax.scatter(x_i[0], x_i[1], x_i[2], color="orange", label="$x_i$")
465
466 ax.scatter(x_j[0], x_j[1], x_j[2], color="red", label="$x_j$")
467
468 draw_arrow(ax, x_i, x_i + force_i * fscale_fact, "green", "$f_i$")
469 draw_arrow(ax, x_i, x_i + force_i_exact * fscale_fact, "red", "$f_i$ (exact)")
470
471 abs_error = np.linalg.norm(force_i - force_i_exact)
472 rel_error = abs_error / np.linalg.norm(force_i_exact)
473
474 draw_aabb(ax, box_A, "blue", 0.1)
475 draw_aabb(ax, box_B, "red", 0.1)
476
477 center_view = (0.5, 0.0, 0.0)
478 view_size = 2.0
479 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
480 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
481 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
482 ax.set_xlabel("X")
483 ax.set_ylabel("Y")
484 ax.set_zlabel("Z")
485
486 return ax, rel_error, abs_error
Now let’s put everything together to get a FMM force We start with the following parameters (see figure below for the representation)
499 s_B = (0, 0, 0)
500 s_A = (1, 0, 0)
501
502 box_B_size = 0.5
503 box_A_size = 0.5
504
505 x_j = (0.2, 0.2, 0.0)
506 x_i = (1.2, 0.2, 0.0)
507 m_j = 1
508 m_i = 1
509
510 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
511 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
512 a_i = (x_i[0] - s_A[0], x_i[1] - s_A[1], x_i[2] - s_A[2])
Q_n_B = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0.008000000000000002, v_001=0.008000000000000002, v_002=0, v_011=0.008000000000000002, v_012=0, v_022=0, v_111=0.008000000000000002, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=0.0016000000000000005, v_0001=0.0016000000000000005, v_0002=0, v_0011=0.0016000000000000005, v_0012=0, v_0022=0, v_0111=0.0016000000000000005, v_0112=0, v_0122=0, v_0222=0, v_1111=0.0016000000000000005, v_1112=0, v_1122=0, v_1222=0, v_2222=0)
)
520 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
521 print("D_n =", D_n)
D_n = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=1, v_1=-0, v_2=-0),
t2=SymTensor3d_2(v_00=2, v_01=-0, v_02=-0, v_11=-1, v_12=0, v_22=-1),
t3=SymTensor3d_3(v_000=6, v_001=0, v_002=0, v_011=-3, v_012=0, v_022=-3, v_111=0, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=24, v_0001=0, v_0002=0, v_0011=-12, v_0012=0, v_0022=-12, v_0111=0, v_0112=0, v_0122=0, v_0222=0, v_1111=9, v_1112=0, v_1122=3, v_1222=0, v_2222=9),
t5=SymTensor3d_5(v_00000=120, v_00001=-0, v_00002=-0, v_00011=-60, v_00012=0, v_00022=-60, v_00111=0, v_00112=0, v_00122=0, v_00222=0, v_01111=45, v_01112=-0, v_01122=15, v_01222=-0, v_02222=45, v_11111=0, v_11112=-0, v_11122=0, v_11222=0, v_12222=-0, v_22222=0)
)
524 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
525 print("dM_k =", dM_k)
dM_k = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=1.431, v_1=-0.3600000000000001, v_2=0),
t2=SymTensor3d_2(v_00=-3.3600000000000003, v_01=1.26, v_02=0, v_11=1.56, v_12=-0, v_22=1.8000000000000003),
t3=SymTensor3d_3(v_000=12, v_001=-4.800000000000001, v_002=0, v_011=-5.7, v_012=0, v_022=-6.300000000000001, v_111=3.6000000000000005, v_112=0, v_122=1.2000000000000002, v_222=0),
t4=SymTensor3d_4(v_0000=-48, v_0001=12, v_0002=-0, v_0011=24, v_0012=-0, v_0022=24, v_0111=-9, v_0112=-0, v_0122=-3, v_0222=-0, v_1111=-18, v_1112=-0, v_1122=-6, v_1222=-0, v_2222=-18),
t5=SymTensor3d_5(v_00000=120, v_00001=-0, v_00002=-0, v_00011=-60, v_00012=0, v_00022=-60, v_00111=0, v_00112=0, v_00122=0, v_00222=0, v_01111=45, v_01112=-0, v_01122=15, v_01222=-0, v_02222=45, v_11111=0, v_11112=-0, v_11122=0, v_11222=0, v_12222=-0, v_22222=0)
)
528 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_i)
529 print("a_k =", a_k)
a_k = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=0.19999999999999996, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.03999999999999998, v_01=0.039999999999999994, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0.007999999999999995, v_001=0.007999999999999997, v_002=0, v_011=0.008, v_012=0, v_022=0, v_111=0.008000000000000002, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=0.0015999999999999986, v_0001=0.001599999999999999, v_0002=0, v_0011=0.0015999999999999996, v_0012=0, v_0022=0, v_0111=0.0016, v_0112=0, v_0122=0, v_0222=0, v_1111=0.0016000000000000005, v_1112=0, v_1122=0, v_1222=0, v_2222=0)
)
force_i = [-1.00000000e+00 1.46584134e-16 -0.00000000e+00]
Now we just need the analytical force to compare
540 def analytic_force_i(x_i, x_j, Gconst):
541 force_i_direct = (x_j[0] - x_i[0], x_j[1] - x_i[1], x_j[2] - x_i[2])
542 force_i_direct /= np.linalg.norm(force_i_direct) ** 3
543 force_i_direct *= m_i
544 return force_i_direct
545
546
547 force_i_exact = analytic_force_i(x_i, x_j, Gconst)
548 print("force_i_exact =", force_i_exact)
force_i_exact = [-1. 0. 0.]
This yields the following case
552 ax, rel_error, abs_error = plot_fmm_case(
553 s_A, box_A_size, x_i, s_B, box_B_size, x_j, force_i, force_i_exact, 0.5
554 )
555
556 plt.title(f"FMM, rel error={rel_error}")
557 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
558 plt.show()
559
560 print("force_i =", force_i)
561 print("force_i_exact =", force_i_exact)
562 print("abs error =", abs_error)
563 print("rel error =", rel_error)
force_i = [-1.00000000e+00 1.46584134e-16 -0.00000000e+00]
force_i_exact = [-1. 0. 0.]
abs error = 1.838827341066391e-16
rel error = 1.838827341066391e-16
And yeah the error is insanelly low, but it is the special case where \(a_i = b_j\). Anyway now let’s wrap all of that mess into a function that does it all and see how the error changes depending on the configure and order of the expansion.
FMM in box#
The following is the function to do the same as above but for whatever order
576 def run_fmm(x_i, x_j, s_A, s_B, m_j, order, do_print):
577
578 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
579 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
580 a_i = (x_i[0] - s_A[0], x_i[1] - s_A[1], x_i[2] - s_A[2])
581
582 if do_print:
583 print("x_i =", x_i)
584 print("x_j =", x_j)
585 print("s_A =", s_A)
586 print("s_B =", s_B)
587 print("b_j =", b_j)
588 print("r =", r)
589 print("a_i =", a_i)
590
591 # compute the tensor product of the displacment
592 if order == 1:
593 Q_n_B = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_j)
594 elif order == 2:
595 Q_n_B = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_j)
596 elif order == 3:
597 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
598 elif order == 4:
599 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
600 elif order == 5:
601 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
602 else:
603 raise ValueError("Invalid order")
604
605 # multiply by mass to get the mass moment
606 Q_n_B *= m_j
607
608 if do_print:
609 print("Q_n_B =", Q_n_B)
610
611 # green function gradients
612 if order == 1:
613 D_n = shamrock.phys.green_func_grav_cartesian_1_1(r)
614 elif order == 2:
615 D_n = shamrock.phys.green_func_grav_cartesian_1_2(r)
616 elif order == 3:
617 D_n = shamrock.phys.green_func_grav_cartesian_1_3(r)
618 elif order == 4:
619 D_n = shamrock.phys.green_func_grav_cartesian_1_4(r)
620 elif order == 5:
621 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
622 else:
623 raise ValueError("Invalid order")
624
625 if do_print:
626 print("D_n =", D_n)
627
628 if order == 1:
629 dM_k = shamrock.phys.get_dM_mat_1(D_n, Q_n_B)
630 elif order == 2:
631 dM_k = shamrock.phys.get_dM_mat_2(D_n, Q_n_B)
632 elif order == 3:
633 dM_k = shamrock.phys.get_dM_mat_3(D_n, Q_n_B)
634 elif order == 4:
635 dM_k = shamrock.phys.get_dM_mat_4(D_n, Q_n_B)
636 elif order == 5:
637 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
638 else:
639 raise ValueError("Invalid order")
640
641 if do_print:
642 print("dM_k =", dM_k)
643
644 if order == 1:
645 a_k = shamrock.math.SymTensorCollection_f64_0_0.from_vec(a_i)
646 elif order == 2:
647 a_k = shamrock.math.SymTensorCollection_f64_0_1.from_vec(a_i)
648 elif order == 3:
649 a_k = shamrock.math.SymTensorCollection_f64_0_2.from_vec(a_i)
650 elif order == 4:
651 a_k = shamrock.math.SymTensorCollection_f64_0_3.from_vec(a_i)
652 elif order == 5:
653 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_i)
654 else:
655 raise ValueError("Invalid order")
656
657 if do_print:
658 print("a_k =", a_k)
659
660 if order == 1:
661 result = shamrock.phys.contract_grav_moment_to_force_1(a_k, dM_k)
662 elif order == 2:
663 result = shamrock.phys.contract_grav_moment_to_force_2(a_k, dM_k)
664 elif order == 3:
665 result = shamrock.phys.contract_grav_moment_to_force_3(a_k, dM_k)
666 elif order == 4:
667 result = shamrock.phys.contract_grav_moment_to_force_4(a_k, dM_k)
668 elif order == 5:
669 result = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
670 else:
671 raise ValueError("Invalid order")
672
673 Gconst = 1 # let's just set the grav constant to 1
674 force_i = -Gconst * np.array(result)
675 if do_print:
676 print("force_i =", force_i)
677
678 force_i_exact = analytic_force_i(x_i, x_j, Gconst)
679 if do_print:
680 print("force_i_exact =", force_i_exact)
681
682 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
683 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
684 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
685
686 angle = (b_A_size + b_B_size) / b_dist
687
688 if do_print:
689 print("b_A_size =", b_A_size)
690 print("b_B_size =", b_B_size)
691 print("b_dist =", b_dist)
692 print("angle =", angle)
693
694 return force_i, force_i_exact, angle
Let’s try with some new parameters
699 s_B = (0, 0, 0)
700 s_A = (1, 0, 0)
701
702 box_B_size = 0.5
703 box_A_size = 0.5
704
705 x_j = (0.2, 0.2, 0.0)
706 x_i = (1.2, 0.2, 0.2)
707 m_j = 1
708 m_i = 1
709
710 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order=5, do_print=True)
711 ax, rel_error, abs_error = plot_fmm_case(
712 s_A, box_A_size, x_i, s_B, box_B_size, x_j, force_i, force_i_exact, 0.5
713 )
714
715 plt.title(f"FMM angle={angle:.5f} rel error={rel_error:.2e}")
716 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
717 plt.show()
718
719 print("force_i =", force_i)
720 print("force_i_exact =", force_i_exact)
721 print("abs error =", abs_error)
722 print("rel error =", rel_error)
x_i = (1.2, 0.2, 0.2)
x_j = (0.2, 0.2, 0.0)
s_A = (1, 0, 0)
s_B = (0, 0, 0)
b_j = (0.2, 0.2, 0.0)
r = (-1, 0, 0)
a_i = (0.19999999999999996, 0.2, 0.2)
Q_n_B = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0.008000000000000002, v_001=0.008000000000000002, v_002=0, v_011=0.008000000000000002, v_012=0, v_022=0, v_111=0.008000000000000002, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=0.0016000000000000005, v_0001=0.0016000000000000005, v_0002=0, v_0011=0.0016000000000000005, v_0012=0, v_0022=0, v_0111=0.0016000000000000005, v_0112=0, v_0122=0, v_0222=0, v_1111=0.0016000000000000005, v_1112=0, v_1122=0, v_1222=0, v_2222=0)
)
D_n = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=1, v_1=-0, v_2=-0),
t2=SymTensor3d_2(v_00=2, v_01=-0, v_02=-0, v_11=-1, v_12=0, v_22=-1),
t3=SymTensor3d_3(v_000=6, v_001=0, v_002=0, v_011=-3, v_012=0, v_022=-3, v_111=0, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=24, v_0001=0, v_0002=0, v_0011=-12, v_0012=0, v_0022=-12, v_0111=0, v_0112=0, v_0122=0, v_0222=0, v_1111=9, v_1112=0, v_1122=3, v_1222=0, v_2222=9),
t5=SymTensor3d_5(v_00000=120, v_00001=-0, v_00002=-0, v_00011=-60, v_00012=0, v_00022=-60, v_00111=0, v_00112=0, v_00122=0, v_00222=0, v_01111=45, v_01112=-0, v_01122=15, v_01222=-0, v_02222=45, v_11111=0, v_11112=-0, v_11122=0, v_11222=0, v_12222=-0, v_22222=0)
)
dM_k = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=1.431, v_1=-0.3600000000000001, v_2=0),
t2=SymTensor3d_2(v_00=-3.3600000000000003, v_01=1.26, v_02=0, v_11=1.56, v_12=-0, v_22=1.8000000000000003),
t3=SymTensor3d_3(v_000=12, v_001=-4.800000000000001, v_002=0, v_011=-5.7, v_012=0, v_022=-6.300000000000001, v_111=3.6000000000000005, v_112=0, v_122=1.2000000000000002, v_222=0),
t4=SymTensor3d_4(v_0000=-48, v_0001=12, v_0002=-0, v_0011=24, v_0012=-0, v_0022=24, v_0111=-9, v_0112=-0, v_0122=-3, v_0222=-0, v_1111=-18, v_1112=-0, v_1122=-6, v_1222=-0, v_2222=-18),
t5=SymTensor3d_5(v_00000=120, v_00001=-0, v_00002=-0, v_00011=-60, v_00012=0, v_00022=-60, v_00111=0, v_00112=0, v_00122=0, v_00222=0, v_01111=45, v_01112=-0, v_01122=15, v_01222=-0, v_02222=45, v_11111=0, v_11112=-0, v_11122=0, v_11222=0, v_12222=-0, v_22222=0)
)
a_k = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=0.19999999999999996, v_1=0.2, v_2=0.2),
t2=SymTensor3d_2(v_00=0.03999999999999998, v_01=0.039999999999999994, v_02=0.039999999999999994, v_11=0.04000000000000001, v_12=0.04000000000000001, v_22=0.04000000000000001),
t3=SymTensor3d_3(v_000=0.007999999999999995, v_001=0.007999999999999997, v_002=0.007999999999999997, v_011=0.008, v_012=0.008, v_022=0.008, v_111=0.008000000000000002, v_112=0.008000000000000002, v_122=0.008000000000000002, v_222=0.008000000000000002),
t4=SymTensor3d_4(v_0000=0.0015999999999999986, v_0001=0.001599999999999999, v_0002=0.001599999999999999, v_0011=0.0015999999999999996, v_0012=0.0015999999999999996, v_0022=0.0015999999999999996, v_0111=0.0016, v_0112=0.0016, v_0122=0.0016, v_0222=0.0016, v_1111=0.0016000000000000005, v_1112=0.0016000000000000005, v_1122=0.0016000000000000005, v_1222=0.0016000000000000005, v_2222=0.0016000000000000005)
)
force_i = [-9.43000000e-01 1.31838984e-16 -1.88000000e-01]
force_i_exact = [-0.94286603 0. -0.18857321]
b_A_size = 0.34641016151377546
b_B_size = 0.28284271247461906
b_dist = 1.0
angle = 0.6292528739883945
force_i = [-9.43000000e-01 1.31838984e-16 -1.88000000e-01]
force_i_exact = [-0.94286603 0. -0.18857321]
abs error = 0.0005886534739763063
rel error = 0.0006121996129353586
Varying the order of the expansion#
729 s_B = (0, 0, 0)
730 s_A = (1, 0, 0)
731
732 box_B_size = 0.5
733 box_A_size = 0.5
734
735 x_j = (0.2, 0.2, 0.0)
736 x_i = (0.8, 0.2, 0.2)
737 m_j = 1
738 m_i = 1
739
740
741 for order in range(1, 6):
742 print("--------------------------------")
743 print(f"Running FMM order = {order}")
744 print("--------------------------------")
745
746 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order, do_print=True)
747 ax, rel_error, abs_error = plot_fmm_case(
748 s_A, box_A_size, x_i, s_B, box_B_size, x_j, force_i, force_i_exact, 0.2
749 )
750
751 plt.title(f"FMM order={order} angle={angle:.5f} rel error={rel_error:.2e}")
752 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
753 plt.show()
754
755 print("force_i =", force_i)
756 print("force_i_exact =", force_i_exact)
757 print("abs error =", abs_error)
758 print("rel error =", rel_error)
--------------------------------
Running FMM order = 1
--------------------------------
x_i = (0.8, 0.2, 0.2)
x_j = (0.2, 0.2, 0.0)
s_A = (1, 0, 0)
s_B = (0, 0, 0)
b_j = (0.2, 0.2, 0.0)
r = (-1, 0, 0)
a_i = (-0.19999999999999996, 0.2, 0.2)
Q_n_B = SymTensorCollection_f64_0_0(t0=1)
D_n = SymTensorCollection_f64_1_1(
t1=SymTensor3d_1(v_0=1, v_1=-0, v_2=-0)
)
dM_k = SymTensorCollection_f64_1_1(
t1=SymTensor3d_1(v_0=1, v_1=-0, v_2=-0)
)
a_k = SymTensorCollection_f64_0_0(t0=1)
force_i = [-1. 0. 0.]
force_i_exact = [-2.37170825 0. -0.79056942]
b_A_size = 0.34641016151377546
b_B_size = 0.28284271247461906
b_dist = 1.0
angle = 0.6292528739883945
force_i = [-1. 0. 0.]
force_i_exact = [-2.37170825 0. -0.79056942]
abs error = 1.5832193498525176
rel error = 0.6332877399410073
--------------------------------
Running FMM order = 2
--------------------------------
x_i = (0.8, 0.2, 0.2)
x_j = (0.2, 0.2, 0.0)
s_A = (1, 0, 0)
s_B = (0, 0, 0)
b_j = (0.2, 0.2, 0.0)
r = (-1, 0, 0)
a_i = (-0.19999999999999996, 0.2, 0.2)
Q_n_B = SymTensorCollection_f64_0_1(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0)
)
D_n = SymTensorCollection_f64_1_2(
t1=SymTensor3d_1(v_0=1, v_1=-0, v_2=-0),
t2=SymTensor3d_2(v_00=2, v_01=-0, v_02=-0, v_11=-1, v_12=0, v_22=-1)
)
dM_k = SymTensorCollection_f64_1_2(
t1=SymTensor3d_1(v_0=1.4, v_1=-0.2, v_2=0),
t2=SymTensor3d_2(v_00=-2, v_01=0, v_02=0, v_11=1, v_12=-0, v_22=1)
)
a_k = SymTensorCollection_f64_0_1(
t0=1,
t1=SymTensor3d_1(v_0=-0.19999999999999996, v_1=0.2, v_2=0.2)
)
force_i = [-1.8 -0. -0.2]
force_i_exact = [-2.37170825 0. -0.79056942]
b_A_size = 0.34641016151377546
b_B_size = 0.28284271247461906
b_dist = 1.0
angle = 0.6292528739883945
force_i = [-1.8 -0. -0.2]
force_i_exact = [-2.37170825 0. -0.79056942]
abs error = 0.8219626217344294
rel error = 0.32878504869377184
--------------------------------
Running FMM order = 3
--------------------------------
x_i = (0.8, 0.2, 0.2)
x_j = (0.2, 0.2, 0.0)
s_A = (1, 0, 0)
s_B = (0, 0, 0)
b_j = (0.2, 0.2, 0.0)
r = (-1, 0, 0)
a_i = (-0.19999999999999996, 0.2, 0.2)
Q_n_B = SymTensorCollection_f64_0_2(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0)
)
D_n = SymTensorCollection_f64_1_3(
t1=SymTensor3d_1(v_0=1, v_1=-0, v_2=-0),
t2=SymTensor3d_2(v_00=2, v_01=-0, v_02=-0, v_11=-1, v_12=0, v_22=-1),
t3=SymTensor3d_3(v_000=6, v_001=0, v_002=0, v_011=-3, v_012=0, v_022=-3, v_111=0, v_112=0, v_122=0, v_222=0)
)
dM_k = SymTensorCollection_f64_1_3(
t1=SymTensor3d_1(v_0=1.46, v_1=-0.32000000000000006, v_2=0),
t2=SymTensor3d_2(v_00=-3.2, v_01=0.6000000000000001, v_02=0, v_11=1.6, v_12=-0, v_22=1.6),
t3=SymTensor3d_3(v_000=6, v_001=0, v_002=0, v_011=-3, v_012=0, v_022=-3, v_111=0, v_112=0, v_122=0, v_222=0)
)
a_k = SymTensorCollection_f64_0_2(
t0=1,
t1=SymTensor3d_1(v_0=-0.19999999999999996, v_1=0.2, v_2=0.2),
t2=SymTensor3d_2(v_00=0.03999999999999998, v_01=-0.039999999999999994, v_02=-0.039999999999999994, v_11=0.04000000000000001, v_12=0.04000000000000001, v_22=0.04000000000000001)
)
force_i = [-2.22000000e+00 4.16333634e-17 -4.40000000e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
b_A_size = 0.34641016151377546
b_B_size = 0.28284271247461906
b_dist = 1.0
angle = 0.6292528739883945
force_i = [-2.22000000e+00 4.16333634e-17 -4.40000000e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
abs error = 0.3819873118341143
rel error = 0.15279492473364575
--------------------------------
Running FMM order = 4
--------------------------------
x_i = (0.8, 0.2, 0.2)
x_j = (0.2, 0.2, 0.0)
s_A = (1, 0, 0)
s_B = (0, 0, 0)
b_j = (0.2, 0.2, 0.0)
r = (-1, 0, 0)
a_i = (-0.19999999999999996, 0.2, 0.2)
Q_n_B = SymTensorCollection_f64_0_3(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0.008000000000000002, v_001=0.008000000000000002, v_002=0, v_011=0.008000000000000002, v_012=0, v_022=0, v_111=0.008000000000000002, v_112=0, v_122=0, v_222=0)
)
D_n = SymTensorCollection_f64_1_4(
t1=SymTensor3d_1(v_0=1, v_1=-0, v_2=-0),
t2=SymTensor3d_2(v_00=2, v_01=-0, v_02=-0, v_11=-1, v_12=0, v_22=-1),
t3=SymTensor3d_3(v_000=6, v_001=0, v_002=0, v_011=-3, v_012=0, v_022=-3, v_111=0, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=24, v_0001=0, v_0002=0, v_0011=-12, v_0012=0, v_0022=-12, v_0111=0, v_0112=0, v_0122=0, v_0222=0, v_1111=9, v_1112=0, v_1122=3, v_1222=0, v_2222=9)
)
dM_k = SymTensorCollection_f64_1_4(
t1=SymTensor3d_1(v_0=1.444, v_1=-0.3560000000000001, v_2=0),
t2=SymTensor3d_2(v_00=-3.4400000000000004, v_01=1.08, v_02=0, v_11=1.6600000000000001, v_12=-0, v_22=1.7800000000000002),
t3=SymTensor3d_3(v_000=10.8, v_001=-2.4000000000000004, v_002=0, v_011=-5.4, v_012=0, v_022=-5.4, v_111=1.8, v_112=0, v_122=0.6000000000000001, v_222=0),
t4=SymTensor3d_4(v_0000=-24, v_0001=-0, v_0002=-0, v_0011=12, v_0012=-0, v_0022=12, v_0111=-0, v_0112=-0, v_0122=-0, v_0222=-0, v_1111=-9, v_1112=-0, v_1122=-3, v_1222=-0, v_2222=-9)
)
a_k = SymTensorCollection_f64_0_3(
t0=1,
t1=SymTensor3d_1(v_0=-0.19999999999999996, v_1=0.2, v_2=0.2),
t2=SymTensor3d_2(v_00=0.03999999999999998, v_01=-0.039999999999999994, v_02=-0.039999999999999994, v_11=0.04000000000000001, v_12=0.04000000000000001, v_22=0.04000000000000001),
t3=SymTensor3d_3(v_000=-0.007999999999999995, v_001=0.007999999999999997, v_002=0.007999999999999997, v_011=-0.008, v_012=-0.008, v_022=-0.008, v_111=0.008000000000000002, v_112=0.008000000000000002, v_122=0.008000000000000002, v_222=0.008000000000000002)
)
force_i = [-2.38000000e+00 4.51028104e-17 -6.20000000e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
b_A_size = 0.34641016151377546
b_B_size = 0.28284271247461906
b_dist = 1.0
angle = 0.6292528739883945
force_i = [-2.38000000e+00 4.51028104e-17 -6.20000000e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
abs error = 0.17077083634710005
rel error = 0.06830833453884004
--------------------------------
Running FMM order = 5
--------------------------------
x_i = (0.8, 0.2, 0.2)
x_j = (0.2, 0.2, 0.0)
s_A = (1, 0, 0)
s_B = (0, 0, 0)
b_j = (0.2, 0.2, 0.0)
r = (-1, 0, 0)
a_i = (-0.19999999999999996, 0.2, 0.2)
Q_n_B = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0.008000000000000002, v_001=0.008000000000000002, v_002=0, v_011=0.008000000000000002, v_012=0, v_022=0, v_111=0.008000000000000002, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=0.0016000000000000005, v_0001=0.0016000000000000005, v_0002=0, v_0011=0.0016000000000000005, v_0012=0, v_0022=0, v_0111=0.0016000000000000005, v_0112=0, v_0122=0, v_0222=0, v_1111=0.0016000000000000005, v_1112=0, v_1122=0, v_1222=0, v_2222=0)
)
D_n = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=1, v_1=-0, v_2=-0),
t2=SymTensor3d_2(v_00=2, v_01=-0, v_02=-0, v_11=-1, v_12=0, v_22=-1),
t3=SymTensor3d_3(v_000=6, v_001=0, v_002=0, v_011=-3, v_012=0, v_022=-3, v_111=0, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=24, v_0001=0, v_0002=0, v_0011=-12, v_0012=0, v_0022=-12, v_0111=0, v_0112=0, v_0122=0, v_0222=0, v_1111=9, v_1112=0, v_1122=3, v_1222=0, v_2222=9),
t5=SymTensor3d_5(v_00000=120, v_00001=-0, v_00002=-0, v_00011=-60, v_00012=0, v_00022=-60, v_00111=0, v_00112=0, v_00122=0, v_00222=0, v_01111=45, v_01112=-0, v_01122=15, v_01222=-0, v_02222=45, v_11111=0, v_11112=-0, v_11122=0, v_11222=0, v_12222=-0, v_22222=0)
)
dM_k = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=1.431, v_1=-0.3600000000000001, v_2=0),
t2=SymTensor3d_2(v_00=-3.3600000000000003, v_01=1.26, v_02=0, v_11=1.56, v_12=-0, v_22=1.8000000000000003),
t3=SymTensor3d_3(v_000=12, v_001=-4.800000000000001, v_002=0, v_011=-5.7, v_012=0, v_022=-6.300000000000001, v_111=3.6000000000000005, v_112=0, v_122=1.2000000000000002, v_222=0),
t4=SymTensor3d_4(v_0000=-48, v_0001=12, v_0002=-0, v_0011=24, v_0012=-0, v_0022=24, v_0111=-9, v_0112=-0, v_0122=-3, v_0222=-0, v_1111=-18, v_1112=-0, v_1122=-6, v_1222=-0, v_2222=-18),
t5=SymTensor3d_5(v_00000=120, v_00001=-0, v_00002=-0, v_00011=-60, v_00012=0, v_00022=-60, v_00111=0, v_00112=0, v_00122=0, v_00222=0, v_01111=45, v_01112=-0, v_01122=15, v_01222=-0, v_02222=45, v_11111=0, v_11112=-0, v_11122=0, v_11222=0, v_12222=-0, v_22222=0)
)
a_k = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=-0.19999999999999996, v_1=0.2, v_2=0.2),
t2=SymTensor3d_2(v_00=0.03999999999999998, v_01=-0.039999999999999994, v_02=-0.039999999999999994, v_11=0.04000000000000001, v_12=0.04000000000000001, v_22=0.04000000000000001),
t3=SymTensor3d_3(v_000=-0.007999999999999995, v_001=0.007999999999999997, v_002=0.007999999999999997, v_011=-0.008, v_012=-0.008, v_022=-0.008, v_111=0.008000000000000002, v_112=0.008000000000000002, v_122=0.008000000000000002, v_222=0.008000000000000002),
t4=SymTensor3d_4(v_0000=0.0015999999999999986, v_0001=-0.001599999999999999, v_0002=-0.001599999999999999, v_0011=0.0015999999999999996, v_0012=0.0015999999999999996, v_0022=0.0015999999999999996, v_0111=-0.0016, v_0112=-0.0016, v_0122=-0.0016, v_0222=-0.0016, v_1111=0.0016000000000000005, v_1112=0.0016000000000000005, v_1122=0.0016000000000000005, v_1222=0.0016000000000000005, v_2222=0.0016000000000000005)
)
force_i = [-2.41500000e+00 3.12250226e-17 -7.24000000e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
b_A_size = 0.34641016151377546
b_B_size = 0.28284271247461906
b_dist = 1.0
angle = 0.6292528739883945
force_i = [-2.41500000e+00 3.12250226e-17 -7.24000000e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
abs error = 0.07940820523782492
rel error = 0.03176328209512998
Sweeping through angles#
766 s_B = (0, 0, 0)
767 s_A_all = [(0.8, 0, 0), (1, 0, 0), (1.2, 0, 0)]
768
769 box_B_size = 0.5
770 box_A_size = 0.5
771
772 x_j = (0.2, 0.2, 0.0)
773 x_i = (0.8, 0.2, 0.2)
774 m_j = 1
775 m_i = 1
776
777 order = 3
778
779 for s_A in s_A_all:
780 print("--------------------------------")
781 print(f"Running FMM s_a = {s_A}")
782 print("--------------------------------")
783
784 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order, do_print=True)
785 ax, rel_error, abs_error = plot_fmm_case(
786 s_A, box_A_size, x_i, s_B, box_B_size, x_j, force_i, force_i_exact, 0.2
787 )
788
789 plt.title(f"FMM order={order} angle={angle:.5f} rel error={rel_error:.2e}")
790 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
791 plt.show()
792
793 print("force_i =", force_i)
794 print("force_i_exact =", force_i_exact)
795 print("abs error =", abs_error)
796 print("rel error =", rel_error)
--------------------------------
Running FMM s_a = (0.8, 0, 0)
--------------------------------
x_i = (0.8, 0.2, 0.2)
x_j = (0.2, 0.2, 0.0)
s_A = (0.8, 0, 0)
s_B = (0, 0, 0)
b_j = (0.2, 0.2, 0.0)
r = (-0.8, 0, 0)
a_i = (0.0, 0.2, 0.2)
Q_n_B = SymTensorCollection_f64_0_2(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0)
)
D_n = SymTensorCollection_f64_1_3(
t1=SymTensor3d_1(v_0=1.5625, v_1=-0, v_2=-0),
t2=SymTensor3d_2(v_00=3.906249999999999, v_01=-0, v_02=-0, v_11=-1.9531249999999998, v_12=0, v_22=-1.9531249999999998),
t3=SymTensor3d_3(v_000=14.648437499999995, v_001=0, v_002=0, v_011=-7.324218749999997, v_012=0, v_022=-7.324218749999997, v_111=0, v_112=0, v_122=0, v_222=0)
)
dM_k = SymTensorCollection_f64_1_3(
t1=SymTensor3d_1(v_0=2.490234375, v_1=-0.68359375, v_2=0),
t2=SymTensor3d_2(v_00=-6.835937499999998, v_01=1.4648437499999996, v_02=0, v_11=3.417968749999999, v_12=-0, v_22=3.417968749999999),
t3=SymTensor3d_3(v_000=14.648437499999995, v_001=0, v_002=0, v_011=-7.324218749999997, v_012=0, v_022=-7.324218749999997, v_111=0, v_112=0, v_122=0, v_222=0)
)
a_k = SymTensorCollection_f64_0_2(
t0=1,
t1=SymTensor3d_1(v_0=0, v_1=0.2, v_2=0.2),
t2=SymTensor3d_2(v_00=0, v_01=0, v_02=0, v_11=0.04000000000000001, v_12=0.04000000000000001, v_22=0.04000000000000001)
)
force_i = [-2.49023438e+00 1.11022302e-16 -6.83593750e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
b_A_size = 0.28284271247461906
b_B_size = 0.28284271247461906
b_dist = 0.8
angle = 0.7071067811865476
force_i = [-2.49023438e+00 1.11022302e-16 -6.83593750e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
abs error = 0.15966288352037078
rel error = 0.06386515340814833
--------------------------------
Running FMM s_a = (1, 0, 0)
--------------------------------
x_i = (0.8, 0.2, 0.2)
x_j = (0.2, 0.2, 0.0)
s_A = (1, 0, 0)
s_B = (0, 0, 0)
b_j = (0.2, 0.2, 0.0)
r = (-1, 0, 0)
a_i = (-0.19999999999999996, 0.2, 0.2)
Q_n_B = SymTensorCollection_f64_0_2(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0)
)
D_n = SymTensorCollection_f64_1_3(
t1=SymTensor3d_1(v_0=1, v_1=-0, v_2=-0),
t2=SymTensor3d_2(v_00=2, v_01=-0, v_02=-0, v_11=-1, v_12=0, v_22=-1),
t3=SymTensor3d_3(v_000=6, v_001=0, v_002=0, v_011=-3, v_012=0, v_022=-3, v_111=0, v_112=0, v_122=0, v_222=0)
)
dM_k = SymTensorCollection_f64_1_3(
t1=SymTensor3d_1(v_0=1.46, v_1=-0.32000000000000006, v_2=0),
t2=SymTensor3d_2(v_00=-3.2, v_01=0.6000000000000001, v_02=0, v_11=1.6, v_12=-0, v_22=1.6),
t3=SymTensor3d_3(v_000=6, v_001=0, v_002=0, v_011=-3, v_012=0, v_022=-3, v_111=0, v_112=0, v_122=0, v_222=0)
)
a_k = SymTensorCollection_f64_0_2(
t0=1,
t1=SymTensor3d_1(v_0=-0.19999999999999996, v_1=0.2, v_2=0.2),
t2=SymTensor3d_2(v_00=0.03999999999999998, v_01=-0.039999999999999994, v_02=-0.039999999999999994, v_11=0.04000000000000001, v_12=0.04000000000000001, v_22=0.04000000000000001)
)
force_i = [-2.22000000e+00 4.16333634e-17 -4.40000000e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
b_A_size = 0.34641016151377546
b_B_size = 0.28284271247461906
b_dist = 1.0
angle = 0.6292528739883945
force_i = [-2.22000000e+00 4.16333634e-17 -4.40000000e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
abs error = 0.3819873118341143
rel error = 0.15279492473364575
--------------------------------
Running FMM s_a = (1.2, 0, 0)
--------------------------------
x_i = (0.8, 0.2, 0.2)
x_j = (0.2, 0.2, 0.0)
s_A = (1.2, 0, 0)
s_B = (0, 0, 0)
b_j = (0.2, 0.2, 0.0)
r = (-1.2, 0, 0)
a_i = (-0.3999999999999999, 0.2, 0.2)
Q_n_B = SymTensorCollection_f64_0_2(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0)
)
D_n = SymTensorCollection_f64_1_3(
t1=SymTensor3d_1(v_0=0.6944444444444444, v_1=-0, v_2=-0),
t2=SymTensor3d_2(v_00=1.157407407407407, v_01=-0, v_02=-0, v_11=-0.5787037037037037, v_12=0, v_22=-0.5787037037037037),
t3=SymTensor3d_3(v_000=2.893518518518519, v_001=0, v_002=0, v_011=-1.4467592592592589, v_012=0, v_022=-1.4467592592592589, v_111=0, v_112=0, v_122=0, v_222=0)
)
dM_k = SymTensorCollection_f64_1_3(
t1=SymTensor3d_1(v_0=0.954861111111111, v_1=-0.1736111111111111, v_2=0),
t2=SymTensor3d_2(v_00=-1.7361111111111107, v_01=0.2893518518518518, v_02=0, v_11=0.8680555555555556, v_12=-0, v_22=0.8680555555555556),
t3=SymTensor3d_3(v_000=2.893518518518519, v_001=0, v_002=0, v_011=-1.4467592592592589, v_012=0, v_022=-1.4467592592592589, v_111=0, v_112=0, v_122=0, v_222=0)
)
a_k = SymTensorCollection_f64_0_2(
t0=1,
t1=SymTensor3d_1(v_0=-0.3999999999999999, v_1=0.2, v_2=0.2),
t2=SymTensor3d_2(v_00=0.15999999999999992, v_01=-0.07999999999999999, v_02=-0.07999999999999999, v_11=0.04000000000000001, v_12=0.04000000000000001, v_22=0.04000000000000001)
)
force_i = [-1.88078704e+00 -1.38777878e-17 -2.89351852e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
b_A_size = 0.48989794855663554
b_B_size = 0.28284271247461906
b_dist = 1.2
angle = 0.6439505508593789
force_i = [-1.88078704e+00 -1.38777878e-17 -2.89351852e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
abs error = 0.7015858309588148
rel error = 0.28063433238352603
FMM precision (Angle)#
For this test we will generate a pair of random positions \(x_i\) and \(x_j\). Then we will generate two boxes around the positions \(s_A\) and \(s_B\) where each is at a distance box_scale_fact from their respective particle. We then perform the FMM expansion to compute the force on \(x_i\) as well as the exact force. We will then plot the relative error as a function of the angle \(\theta = (b_A + b_B) / |\mathbf{s}_A - \mathbf{s}_B|\) where \(b_A\) and \(b_B\) are the distances from the particle to the box centers.
810 plt.figure()
811 for order in range(1, 6):
812 print("--------------------------------")
813 print(f"Running FMM order = {order}")
814 print("--------------------------------")
815
816 # set seed
817 rng = np.random.default_rng(111)
818
819 N = 50000
820
821 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
822 x_i_all = []
823 for i in range(N):
824 x_i_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
825
826 # same for x_j
827 x_j_all = []
828 for i in range(N):
829 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
830
831 box_scale_fact_all = np.linspace(0, 0.1, N).tolist()
832
833 # same for box_1_center
834 s_A_all = []
835 for p, box_scale_fact in zip(x_i_all, box_scale_fact_all):
836 s_A_all.append(
837 (
838 p[0] + box_scale_fact * rng.uniform(-1, 1),
839 p[1] + box_scale_fact * rng.uniform(-1, 1),
840 p[2] + box_scale_fact * rng.uniform(-1, 1),
841 )
842 )
843
844 # same for box_2_center
845 s_B_all = []
846 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
847 s_B_all.append(
848 (
849 p[0] + box_scale_fact * rng.uniform(-1, 1),
850 p[1] + box_scale_fact * rng.uniform(-1, 1),
851 p[2] + box_scale_fact * rng.uniform(-1, 1),
852 )
853 )
854
855 angles = []
856 rel_errors = []
857
858 for x_i, x_j, s_A, s_B in zip(x_i_all, x_j_all, s_A_all, s_B_all):
859
860 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order, do_print=False)
861
862 abs_error = np.linalg.norm(force_i - force_i_exact)
863 rel_error = abs_error / np.linalg.norm(force_i_exact)
864
865 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
866 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
867 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
868 angle = (b_A_size + b_B_size) / b_dist
869
870 if angle > 5.0 or angle < 1e-4:
871 continue
872
873 angles.append(angle)
874 rel_errors.append(rel_error)
875
876 print(f"Computed for {len(angles)} cases")
877
878 plt.scatter(angles, rel_errors, s=1, label=f"FMM order = {order}")
879
880
881 def plot_powerlaw(order, center_y):
882 X = [1e-3, 1e-2 / 3, 1e-1]
883 Y = [center_y * (x) ** order for x in X]
884 plt.plot(X, Y, linestyle="dashed", color="black")
885 bbox = dict(boxstyle="round", fc="blanchedalmond", ec="orange", alpha=0.9)
886 plt.text(X[1], Y[1], f"$\\propto x^{order}$", fontsize=9, bbox=bbox)
887
888
889 plot_powerlaw(1, 1)
890 plot_powerlaw(2, 1)
891 plot_powerlaw(3, 1)
892 plot_powerlaw(4, 1)
893 plot_powerlaw(5, 1)
894
895 plt.xlabel("Angle")
896 plt.ylabel("Relative Error")
897 plt.xscale("log")
898 plt.yscale("log")
899 plt.title("FMM precision")
900 plt.legend(loc="lower right")
901 plt.grid()
902 plt.show()
--------------------------------
Running FMM order = 1
--------------------------------
Computed for 49962 cases
--------------------------------
Running FMM order = 2
--------------------------------
Computed for 49962 cases
--------------------------------
Running FMM order = 3
--------------------------------
Computed for 49962 cases
--------------------------------
Running FMM order = 4
--------------------------------
Computed for 49962 cases
--------------------------------
Running FMM order = 5
--------------------------------
Computed for 49962 cases
Mass moment offset#
Now that we know how to compute a FMM force, we now need some remaining tools to exploit it fully in a code. In a code using a tree the procedure to using a FMM is to first propagate the mass moment upward from leafs cells up to the root. Then compute the gravitation moments for all cell-cell interations and then propagate the gravitational moment downward down to the leaves.
We start with the upward pass for the mass moment. To perform it we need to compute the mass moment of a parent according to the one of its children. The issue is that the childrens and the parents do not share the same center. Therefor we need to offset the mass moment of the children to the parent center before summing their moments to get the parent’s one.
This is what we call mass moment translation/offset. This section will showcase its usage and precision.
We start of by defining a particle \(x_j\) and a box \(s_B\) around it as well as a new box \(s_B'\). The goal will be to offset the mass moment of the box \(s_B\) to the box \(s_B'\) and compare it to the moment of the box \(s_B'\) computed directly. This should yield the same result meaning that we never need to compute the moment directly at the parent center and simply use its childrens instead.
932 s_B = (0, 0, 0)
933 box_B_size = 1.0
934 x_j = (0.2, 0.2, 0.0)
935 m_j = 1
936
937 s_B_new = (0.3, 0.3, 0.3)
def plot_mass_moment_offset(s_B, s_B_new, box_B_size):
946 def plot_mass_moment_offset(s_B, s_B_new, box_B_size):
947 box_B = shamrock.math.AABB_f64_3(
948 (
949 s_B[0] - box_B_size / 2,
950 s_B[1] - box_B_size / 2,
951 s_B[2] - box_B_size / 2,
952 ),
953 (
954 s_B[0] + box_B_size / 2,
955 s_B[1] + box_B_size / 2,
956 s_B[2] + box_B_size / 2,
957 ),
958 )
959
960 box_B_new = shamrock.math.AABB_f64_3(
961 (
962 s_B_new[0] - box_B_size / 2,
963 s_B_new[1] - box_B_size / 2,
964 s_B_new[2] - box_B_size / 2,
965 ),
966 (
967 s_B_new[0] + box_B_size / 2,
968 s_B_new[1] + box_B_size / 2,
969 s_B_new[2] + box_B_size / 2,
970 ),
971 )
972
973 fig = plt.figure()
974 ax = fig.add_subplot(111, projection="3d")
975
976 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
977 draw_arrow(ax, s_B_new, x_j, "red", "$b_j' = x_j - s_B'$")
978
979 ax.scatter(s_B[0], s_B[1], s_B[2], color="black", label="s_B")
980 ax.scatter(s_B_new[0], s_B_new[1], s_B_new[2], color="red", label="s_B'")
981 ax.scatter(x_j[0], x_j[1], x_j[2], color="blue", label="$x_j$")
982
983 draw_aabb(ax, box_B, "blue", 0.2)
984 draw_aabb(ax, box_B_new, "red", 0.2)
985
986 center_view = (0.0, 0.0, 0.0)
987 view_size = 2.0
988 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
989 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
990 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
991 ax.set_xlabel("X")
992 ax.set_ylabel("Y")
993 ax.set_zlabel("Z")
994
995 return ax
1005 plot_mass_moment_offset(s_B, s_B_new, box_B_size)
1006
1007 plt.title("Mass moment offset illustration")
1008 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
1009 plt.show()
Moment for box B
b_j = (0.2, 0.2, 0.0)
Q_n_B = SymTensorCollection_f64_0_5(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0.008000000000000002, v_001=0.008000000000000002, v_002=0, v_011=0.008000000000000002, v_012=0, v_022=0, v_111=0.008000000000000002, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=0.0016000000000000005, v_0001=0.0016000000000000005, v_0002=0, v_0011=0.0016000000000000005, v_0012=0, v_0022=0, v_0111=0.0016000000000000005, v_0112=0, v_0122=0, v_0222=0, v_1111=0.0016000000000000005, v_1112=0, v_1122=0, v_1222=0, v_2222=0),
t5=SymTensor3d_5(v_00000=0.00032000000000000013, v_00001=0.00032000000000000013, v_00002=0, v_00011=0.00032000000000000013, v_00012=0, v_00022=0, v_00111=0.00032000000000000013, v_00112=0, v_00122=0, v_00222=0, v_01111=0.00032000000000000013, v_01112=0, v_01122=0, v_01222=0, v_02222=0, v_11111=0.00032000000000000013, v_11112=0, v_11122=0, v_11222=0, v_12222=0, v_22222=0)
)
Moment for box B’
b_jp = (-0.09999999999999998, -0.09999999999999998, -0.3)
Q_n_Bp = SymTensorCollection_f64_0_5(
t0=1,
t1=SymTensor3d_1(v_0=-0.09999999999999998, v_1=-0.09999999999999998, v_2=-0.3),
t2=SymTensor3d_2(v_00=0.009999999999999995, v_01=0.009999999999999995, v_02=0.029999999999999992, v_11=0.009999999999999995, v_12=0.029999999999999992, v_22=0.09),
t3=SymTensor3d_3(v_000=-0.0009999999999999994, v_001=-0.0009999999999999994, v_002=-0.0029999999999999983, v_011=-0.0009999999999999994, v_012=-0.0029999999999999983, v_022=-0.008999999999999998, v_111=-0.0009999999999999994, v_112=-0.0029999999999999983, v_122=-0.008999999999999998, v_222=-0.027),
t4=SymTensor3d_4(v_0000=9.999999999999991e-05, v_0001=9.999999999999991e-05, v_0002=0.00029999999999999976, v_0011=9.999999999999991e-05, v_0012=0.00029999999999999976, v_0022=0.0008999999999999995, v_0111=9.999999999999991e-05, v_0112=0.00029999999999999976, v_0122=0.0008999999999999995, v_0222=0.0026999999999999993, v_1111=9.999999999999991e-05, v_1112=0.00029999999999999976, v_1122=0.0008999999999999995, v_1222=0.0026999999999999993, v_2222=0.0081),
t5=SymTensor3d_5(v_00000=-9.999999999999989e-06, v_00001=-9.999999999999989e-06, v_00002=-2.999999999999997e-05, v_00011=-9.999999999999989e-06, v_00012=-2.999999999999997e-05, v_00022=-8.999999999999994e-05, v_00111=-9.999999999999989e-06, v_00112=-2.999999999999997e-05, v_00122=-8.999999999999994e-05, v_00222=-0.00026999999999999984, v_01111=-9.999999999999989e-06, v_01112=-2.999999999999997e-05, v_01122=-8.999999999999994e-05, v_01222=-0.00026999999999999984, v_02222=-0.0008099999999999997, v_11111=-9.999999999999989e-06, v_11112=-2.999999999999997e-05, v_11122=-8.999999999999994e-05, v_11222=-0.00026999999999999984, v_12222=-0.0008099999999999997, v_22222=-0.00243)
)
Offset the moment in box B to box B’
1029 Q_n_B_offset = shamrock.phys.offset_multipole_5(Q_n_B, s_B, s_B_new)
1030 print("Q_n_B_offset =", Q_n_B_offset)
Q_n_B_offset = SymTensorCollection_f64_0_5(
t0=1,
t1=SymTensor3d_1(v_0=-0.09999999999999998, v_1=-0.09999999999999998, v_2=-0.3),
t2=SymTensor3d_2(v_00=0.009999999999999998, v_01=0.009999999999999998, v_02=0.029999999999999995, v_11=0.009999999999999998, v_12=0.029999999999999995, v_22=0.09),
t3=SymTensor3d_3(v_000=-0.0010000000000000009, v_001=-0.0010000000000000009, v_002=-0.003000000000000001, v_011=-0.0010000000000000009, v_012=-0.003000000000000001, v_022=-0.009, v_111=-0.0010000000000000009, v_112=-0.003000000000000001, v_122=-0.009000000000000001, v_222=-0.027),
t4=SymTensor3d_4(v_0000=0.00010000000000000048, v_0001=0.00010000000000000048, v_0002=0.0003000000000000008, v_0011=0.00010000000000000048, v_0012=0.0003000000000000012, v_0022=0.0009000000000000006, v_0111=0.00010000000000000048, v_0112=0.0003000000000000012, v_0122=0.0009000000000000012, v_0222=0.0027, v_1111=0.00010000000000000048, v_1112=0.0003000000000000012, v_1122=0.0009000000000000012, v_1222=0.0027, v_2222=0.0081),
t5=SymTensor3d_5(v_00000=-1.0000000000000189e-05, v_00001=-1.0000000000000189e-05, v_00002=-3.0000000000000458e-05, v_00011=-1.0000000000000189e-05, v_00012=-3.0000000000000675e-05, v_00022=-9.000000000000056e-05, v_00111=-1.0000000000000189e-05, v_00112=-3.0000000000000675e-05, v_00122=-9.000000000000067e-05, v_00222=-0.00027000000000000055, v_01111=-1.0000000000000189e-05, v_01112=-3.0000000000000675e-05, v_01122=-9.000000000000067e-05, v_01222=-0.00027000000000000044, v_02222=-0.0008100000000000001, v_11111=-1.0000000000000189e-05, v_11112=-3.0000000000000675e-05, v_11122=-9.000000000000067e-05, v_11222=-0.00027000000000000044, v_12222=-0.0008100000000000001, v_22222=-0.00243)
)
Print the norm of the moment in box B’
1036 def tensor_collect_norm(d):
1037 # detect the type of the tensor collection
1038 if isinstance(d, shamrock.math.SymTensorCollection_f64_0_5):
1039 return (
1040 np.sqrt(d.t0 * d.t0)
1041 + np.sqrt(d.t1.inner(d.t1))
1042 + np.sqrt(d.t2.inner(d.t2)) / 2
1043 + np.sqrt(d.t3.inner(d.t3)) / 6
1044 + np.sqrt(d.t4.inner(d.t4)) / 24
1045 + np.sqrt(d.t5.inner(d.t5)) / 120
1046 )
1047 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_4):
1048 return (
1049 np.sqrt(d.t0 * d.t0)
1050 + np.sqrt(d.t1.inner(d.t1))
1051 + np.sqrt(d.t2.inner(d.t2)) / 2
1052 + np.sqrt(d.t3.inner(d.t3)) / 6
1053 + np.sqrt(d.t4.inner(d.t4)) / 24
1054 )
1055 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_3):
1056 return (
1057 np.sqrt(d.t0 * d.t0)
1058 + np.sqrt(d.t1.inner(d.t1))
1059 + np.sqrt(d.t2.inner(d.t2)) / 2
1060 + np.sqrt(d.t3.inner(d.t3)) / 6
1061 )
1062 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_2):
1063 return np.sqrt(d.t0 * d.t0) + np.sqrt(d.t1.inner(d.t1)) + np.sqrt(d.t2.inner(d.t2)) / 2
1064 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_1):
1065 return np.sqrt(d.t0 * d.t0) + np.sqrt(d.t1.inner(d.t1))
1066 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_0):
1067 return np.sqrt(d.t0 * d.t0)
1068 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_5):
1069 return (
1070 np.sqrt(d.t1.inner(d.t1))
1071 + np.sqrt(d.t2.inner(d.t2)) / 2
1072 + np.sqrt(d.t3.inner(d.t3)) / 6
1073 + np.sqrt(d.t4.inner(d.t4)) / 24
1074 + np.sqrt(d.t5.inner(d.t5)) / 120
1075 )
1076 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_4):
1077 return (
1078 np.sqrt(d.t1.inner(d.t1))
1079 + np.sqrt(d.t2.inner(d.t2)) / 2
1080 + np.sqrt(d.t3.inner(d.t3)) / 6
1081 + np.sqrt(d.t4.inner(d.t4)) / 24
1082 )
1083 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_3):
1084 return (
1085 np.sqrt(d.t1.inner(d.t1))
1086 + np.sqrt(d.t2.inner(d.t2)) / 2
1087 + np.sqrt(d.t3.inner(d.t3)) / 6
1088 )
1089 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_2):
1090 return np.sqrt(d.t1.inner(d.t1)) + np.sqrt(d.t2.inner(d.t2)) / 2
1091 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_1):
1092 return np.sqrt(d.t1.inner(d.t1))
1093 else:
1094 raise ValueError(f"Unsupported tensor collection type: {type(d)}")
1095
1096
1097 print("Q_n_B norm =", tensor_collect_norm(Q_n_B))
1098 print("Q_n_Bp norm =", tensor_collect_norm(Q_n_Bp))
Q_n_B norm = 1.3268957002522794
Q_n_Bp norm = 1.3932805671178279
Compute the delta between the moments
1102 delta = Q_n_B_offset - Q_n_Bp
1103 print("delta =", delta)
1104
1105
1106 sqdist_t0 = delta.t0 * delta.t0
1107 sqdist_t1 = delta.t1.inner(delta.t1)
1108 sqdist_t2 = delta.t2.inner(delta.t2)
1109 sqdist_t3 = delta.t3.inner(delta.t3)
1110 sqdist_t4 = delta.t4.inner(delta.t4)
1111 sqdist_t5 = delta.t5.inner(delta.t5)
1112 print("sqdist_t0 =", sqdist_t0)
1113 print("sqdist_t1 =", sqdist_t1)
1114 print("sqdist_t2 =", sqdist_t2)
1115 print("sqdist_t3 =", sqdist_t3)
1116 print("sqdist_t4 =", sqdist_t4)
1117 print("sqdist_t5 =", sqdist_t5)
1118
1119 norm_delta = (
1120 np.sqrt(sqdist_t0)
1121 + np.sqrt(sqdist_t1)
1122 + np.sqrt(sqdist_t2) / 2
1123 + np.sqrt(sqdist_t3) / 6
1124 + np.sqrt(sqdist_t4) / 24
1125 + np.sqrt(sqdist_t5) / 120
1126 )
1127 print("norm_delta =", norm_delta)
1128
1129 print("rel error =", tensor_collect_norm(delta) / tensor_collect_norm(Q_n_Bp))
delta = SymTensorCollection_f64_0_5(
t0=0,
t1=SymTensor3d_1(v_0=0, v_1=0, v_2=0),
t2=SymTensor3d_2(v_00=3.469446951953614e-18, v_01=3.469446951953614e-18, v_02=3.469446951953614e-18, v_11=3.469446951953614e-18, v_12=3.469446951953614e-18, v_22=0),
t3=SymTensor3d_3(v_000=-1.5178830414797062e-18, v_001=-1.5178830414797062e-18, v_002=-2.6020852139652106e-18, v_011=-1.5178830414797062e-18, v_012=-2.6020852139652106e-18, v_022=-1.734723475976807e-18, v_111=-1.5178830414797062e-18, v_112=-2.6020852139652106e-18, v_122=-3.469446951953614e-18, v_222=0),
t4=SymTensor3d_4(v_0000=5.692061405548898e-19, v_0001=5.692061405548898e-19, v_0002=1.0299920638612292e-18, v_0011=5.692061405548898e-19, v_0012=1.463672932855431e-18, v_0022=1.0842021724855044e-18, v_0111=5.692061405548898e-19, v_0112=1.463672932855431e-18, v_0122=1.6263032587282567e-18, v_0222=8.673617379884035e-19, v_1111=5.692061405548898e-19, v_1112=1.463672932855431e-18, v_1122=1.6263032587282567e-18, v_1222=8.673617379884035e-19, v_2222=0),
t5=SymTensor3d_5(v_00000=-1.9989977555201488e-19, v_00001=-1.9989977555201488e-19, v_00002=-4.87890977618477e-19, v_00011=-1.9989977555201488e-19, v_00012=-7.047314121155779e-19, v_00022=-6.2341624917916505e-19, v_00111=-1.9989977555201488e-19, v_00112=-7.047314121155779e-19, v_00122=-7.318364664277155e-19, v_00222=-7.047314121155779e-19, v_01111=-1.9989977555201488e-19, v_01112=-7.047314121155779e-19, v_01122=-7.318364664277155e-19, v_01222=-5.963111948670274e-19, v_02222=-3.2526065174565133e-19, v_11111=-1.9989977555201488e-19, v_11112=-7.047314121155779e-19, v_11122=-7.318364664277155e-19, v_11222=-5.963111948670274e-19, v_12222=-3.2526065174565133e-19, v_22222=0)
)
sqdist_t0 = 0.0
sqdist_t1 = 0.0
sqdist_t2 = 9.62964972193618e-35
sqdist_t3 = 1.4482090402130582e-34
sqdist_t4 = 1.3009195980550256e-34
sqdist_t5 = 9.778680365101367e-35
norm_delta = 7.469878656304813e-18
rel error = 5.3613599676891896e-18
We now want to explore the precision of the offset as a function of the order & distance
1134 plt.figure()
1135
1136 for order in range(2, 6):
1137 # set seed
1138 rng = np.random.default_rng(111)
1139
1140 N = 50000
1141
1142 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
1143 x_j_all = []
1144 for i in range(N):
1145 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1146
1147 box_scale_fact_all = np.linspace(0, 1, N).tolist()
1148
1149 # same for box_1_center
1150 s_B_all = []
1151 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1152 s_B_all.append(
1153 (
1154 p[0] + box_scale_fact * rng.uniform(-1, 1),
1155 p[1] + box_scale_fact * rng.uniform(-1, 1),
1156 p[2] + box_scale_fact * rng.uniform(-1, 1),
1157 )
1158 )
1159
1160 # same for box_2_center
1161 s_Bp_all = []
1162 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1163 s_Bp_all.append(
1164 (
1165 p[0] + box_scale_fact * rng.uniform(-1, 1),
1166 p[1] + box_scale_fact * rng.uniform(-1, 1),
1167 p[2] + box_scale_fact * rng.uniform(-1, 1),
1168 )
1169 )
1170
1171 center_distances = []
1172 rel_errors = []
1173 for x_j, s_B, s_Bp in zip(x_j_all, s_B_all, s_Bp_all):
1174
1175 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1176 b_jp = (x_j[0] - s_Bp[0], x_j[1] - s_Bp[1], x_j[2] - s_Bp[2])
1177
1178 if order == 5:
1179 Q_n_B = shamrock.math.SymTensorCollection_f64_0_5.from_vec(b_j)
1180 Q_n_B *= m_j
1181
1182 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_5.from_vec(b_jp)
1183 Q_n_Bp *= m_j
1184
1185 Q_n_B_offset = shamrock.phys.offset_multipole_5(Q_n_B, s_B, s_Bp)
1186 elif order == 4:
1187 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
1188 Q_n_B *= m_j
1189
1190 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_jp)
1191 Q_n_Bp *= m_j
1192
1193 Q_n_B_offset = shamrock.phys.offset_multipole_4(Q_n_B, s_B, s_Bp)
1194 elif order == 3:
1195 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
1196 Q_n_B *= m_j
1197
1198 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_jp)
1199 Q_n_Bp *= m_j
1200
1201 Q_n_B_offset = shamrock.phys.offset_multipole_3(Q_n_B, s_B, s_Bp)
1202 elif order == 2:
1203 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
1204 Q_n_B *= m_j
1205
1206 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_jp)
1207 Q_n_Bp *= m_j
1208
1209 Q_n_B_offset = shamrock.phys.offset_multipole_2(Q_n_B, s_B, s_Bp)
1210 else:
1211 raise ValueError(f"Unsupported offset order: {order}")
1212
1213 delta = Q_n_B_offset - Q_n_Bp
1214
1215 rel_error = tensor_collect_norm(delta) / tensor_collect_norm(Q_n_B)
1216 rel_errors.append(rel_error)
1217
1218 center_distances.append(np.linalg.norm(np.array(s_B) - np.array(s_Bp)))
1219
1220 plt.scatter(center_distances, rel_errors, s=1, label=f"multipole order = {order}")
1221
1222 plt.xlabel("$\\vert \\vert s_B - s_B'\\vert \\vert$")
1223 plt.ylabel(
1224 "$\\vert \\vert Q_n(s_B) - Q_n(s_B') \\vert \\vert / \\vert \\vert Q_n(s_B) \\vert \\vert$ (relative error) "
1225 )
1226 plt.xscale("log")
1227 plt.yscale("log")
1228 plt.title("Mass moment offset precision")
1229 plt.legend(loc="lower right")
1230 plt.grid()
1231 plt.show()
As shown the precision is basically the floating point precision. Also as a result we can observe a small precision loss for high orders.
Gravitational moment offset#
Now that we know how to offset the mass moment, we need to offset the gravitational moment. This is required as we will compute gravitational moments for cell-cell interactions, but we still need to propagate that moment from a parent cell to its children until each leaves contains the complete gravitational moment which will be used to compute the resulting force.
We devise a similar setup to the mass moment offset. We define a particle \(x_j\) and a box of center \(s_B\) around it. We then define a box of center \(s_A\) around the particle \(x_i\) as well as a new box of center \(s_A'\).
The goal will be to offset the gravitational moment of the box \(s_A\) to the box \(s_A'\) and then compute the resulting FMM force on \(x_i\) in the new box and compare it to the force given the FMM in the box \(s_A\). If everything is working correctly they should be equals.
1260 def plot_grav_moment_offset(s_A, s_Ap, s_B, box_A_size, box_B_size, x_j):
1261 box_A = shamrock.math.AABB_f64_3(
1262 (
1263 s_A[0] - box_A_size / 2,
1264 s_A[1] - box_A_size / 2,
1265 s_A[2] - box_A_size / 2,
1266 ),
1267 (
1268 s_A[0] + box_A_size / 2,
1269 s_A[1] + box_A_size / 2,
1270 s_A[2] + box_A_size / 2,
1271 ),
1272 )
1273
1274 box_Ap = shamrock.math.AABB_f64_3(
1275 (
1276 s_Ap[0] - box_A_size / 2,
1277 s_Ap[1] - box_A_size / 2,
1278 s_Ap[2] - box_A_size / 2,
1279 ),
1280 (
1281 s_Ap[0] + box_A_size / 2,
1282 s_Ap[1] + box_A_size / 2,
1283 s_Ap[2] + box_A_size / 2,
1284 ),
1285 )
1286
1287 box_B = shamrock.math.AABB_f64_3(
1288 (
1289 s_B[0] - box_B_size / 2,
1290 s_B[1] - box_B_size / 2,
1291 s_B[2] - box_B_size / 2,
1292 ),
1293 (
1294 s_B[0] + box_B_size / 2,
1295 s_B[1] + box_B_size / 2,
1296 s_B[2] + box_B_size / 2,
1297 ),
1298 )
1299
1300 fig = plt.figure()
1301 ax = fig.add_subplot(111, projection="3d")
1302
1303 draw_arrow(ax, s_A, s_B, "purple", "$r = s_B - s_A$")
1304 draw_arrow(ax, s_Ap, s_B, "purple", "$r' = s_B - s_A'$")
1305
1306 ax.scatter(s_A[0], s_A[1], s_A[2], color="black", label="s_A")
1307 ax.scatter(s_Ap[0], s_Ap[1], s_Ap[2], color="black", label="s_Ap")
1308 ax.scatter(s_B[0], s_B[1], s_B[2], color="black", label="s_B")
1309
1310 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
1311
1312 ax.scatter(x_j[0], x_j[1], x_j[2], color="red", label="$x_j$")
1313
1314 draw_aabb(ax, box_A, "blue", 0.1)
1315 draw_aabb(ax, box_Ap, "cyan", 0.1)
1316 draw_aabb(ax, box_B, "red", 0.1)
1317
1318 center_view = (0.5, 0.0, 0.0)
1319 view_size = 2.0
1320 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
1321 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
1322 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
1323 ax.set_xlabel("X")
1324 ax.set_ylabel("Y")
1325 ax.set_zlabel("Z")
1326
1327 return ax, rel_error, abs_error
1328
1329
1330 s_B = (0, -0.2, -0.2)
1331 s_A = (1, 0, 0)
1332 s_Ap = (1.1, 0.1, 0.0)
1333
1334 box_B_size = 0.5
1335 box_A_size = 0.5
1336
1337 x_j = (0.2, 0.0, -0.5)
1338 x_i = (1.2, 0.2, 0.0)
1339 m_j = 1
1340 m_i = 1
1341
1342 plot_grav_moment_offset(s_A, s_Ap, s_B, box_A_size, box_B_size, x_j)
1343
1344 plt.title("Mass moment offset illustration")
1345 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
1346 plt.show()
1347
1348 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1349 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
1350 rp = (s_B[0] - s_Ap[0], s_B[1] - s_Ap[1], s_B[2] - s_Ap[2])
Q_n_B = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=-0.3),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=-0.06, v_11=0.04000000000000001, v_12=-0.06, v_22=0.09),
t3=SymTensor3d_3(v_000=0.008000000000000002, v_001=0.008000000000000002, v_002=-0.012, v_011=0.008000000000000002, v_012=-0.012, v_022=0.018, v_111=0.008000000000000002, v_112=-0.012, v_122=0.018, v_222=-0.027),
t4=SymTensor3d_4(v_0000=0.0016000000000000005, v_0001=0.0016000000000000005, v_0002=-0.0024000000000000002, v_0011=0.0016000000000000005, v_0012=-0.0024000000000000002, v_0022=0.0036, v_0111=0.0016000000000000005, v_0112=-0.0024000000000000002, v_0122=0.0036, v_0222=-0.0054, v_1111=0.0016000000000000005, v_1112=-0.0024000000000000002, v_1122=0.0036, v_1222=-0.0054, v_2222=0.0081)
)
1359 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1360 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
1361 # print("D_n =",D_n)
1362 print("dM_k =", dM_k)
dM_k = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=0.9330692614896893, v_1=-0.00676020905615491, v_2=0.5975337326753848),
t2=SymTensor3d_2(v_00=-1.298112911867752, v_01=0.11610754537124784, v_02=-1.8134165309956711, v_11=1.2391425006660395, v_12=0.007638654300740003, v_22=0.05897041120171296),
t3=SymTensor3d_3(v_000=3.4412986364200324, v_001=-0.2809327303938822, v_002=7.571603890766815, v_011=-4.260501873206062, v_012=0.7061511531350734, v_022=0.819203236786024, v_111=-0.23000836838894867, v_112=-2.368831572596142, v_122=0.5109410987828296, v_222=-5.202772318170674),
t4=SymTensor3d_4(v_0000=-18.407459386049837, v_0001=-6.722015784651176, v_0002=-26.667390903250002, v_0011=15.484401006966687, v_0012=-6.501343549296469, v_0022=2.923058379083166, v_0111=7.061511531350736, v_0112=6.874788870665981, v_0122=-0.33949574669955473, v_0222=19.79260203258403, v_1111=-12.975527438856975, v_1112=3.5477305530103456, v_1122=-2.5088735681097085, v_1222=2.9536129962861244, v_2222=-0.41418481097345605),
t5=SymTensor3d_5(v_00000=48.09523078243683, v_00001=41.53164634624548, v_00002=41.53164634624548, v_00011=-24.04761539121845, v_00012=15.843134845979208, v_00022=-24.04761539121845, v_00111=-33.327165801006274, v_00112=-8.204480545239235, v_00122=-8.204480545239235, v_00222=-33.327165801006274, v_01111=17.54061357947698, v_01112=-7.921567422989607, v_01122=6.507001811741461, v_01222=-7.921567422989607, v_02222=17.54061357947698, v_11111=28.630807971662435, v_11112=3.508122715895396, v_11122=4.696357829343838, v_11222=4.696357829343838, v_12222=3.508122715895396, v_22222=28.630807971662435)
)
1365 Dp_n = shamrock.phys.green_func_grav_cartesian_1_5(rp)
1366 dMp_k = shamrock.phys.get_dM_mat_5(Dp_n, Q_n_B)
1367 # print("Dp_n =",Dp_n)
1368 print("dMp_k =", dMp_k)
dMp_k = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=0.8051957104860917, v_1=0.0859156812420047, v_2=0.45554150189781367),
t2=SymTensor3d_2(v_00=-1.1506350288987228, v_01=-0.18473251757367656, v_02=-1.2544427285969801, v_11=0.9166504977909566, v_12=-0.14449625923034512, v_22=0.2339845311077664),
t3=SymTensor3d_3(v_000=2.785021505417446, v_001=0.8733332982608633, v_002=4.485575659227913, v_011=-2.6188090665288035, v_012=0.9963229064995676, v_022=-0.166212438888643, v_111=-1.0489286383915493, v_112=-1.2271271607572491, v_122=0.1755953401306858, v_222=-3.2584484984706643),
t4=SymTensor3d_4(v_0000=-10.165140797758387, v_0001=-7.280740259644215, v_0002=-13.378828262430224, v_0011=7.011559807338903, v_0012=-5.056011305707431, v_0022=3.153580990419484, v_0111=6.6880614626113, v_0112=2.736559558321528, v_0122=0.592678797032914, v_0222=10.642268704108691, v_1111=-5.468435464111547, v_1112=2.59084325799438, v_1122=-1.5431243432273574, v_1222=2.4651680477130498, v_2222=-1.6104566471921278),
t5=SymTensor3d_5(v_00000=18.72060355756975, v_00001=26.56716119514329, v_00002=17.71144079676219, v_00011=-5.393848537488466, v_00012=9.5194877791114, v_00022=-13.3267550200813, v_00111=-20.376497872789148, v_00112=-2.1221000974285835, v_00122=-6.190663322354132, v_00222=-15.5893406993336, v_01111=2.744597589049272, v_01112=-4.447970272335893, v_01122=2.649250948439194, v_01222=-5.071517506775504, v_02222=10.677504071642106, v_11111=17.09788111542832, v_11112=0.499017743463504, v_11122=3.2786167573608305, v_11222=1.6230823539650792, v_12222=2.9120465649933016, v_22222=13.966258345368523)
)
Offset the grav moment to dMp_k
1372 dM_k_offset = shamrock.phys.offset_dM_mat_5(dM_k, s_A, s_Ap)
1373 print("dM_k_offset =", dM_k_offset)
dM_k_offset = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=0.8102626162966926, v_1=0.09094581810461017, v_2=0.44647934387825156),
t2=SymTensor3d_2(v_00=-1.0527084113659793, v_01=-0.20033644012740712, v_02=-1.140196465290454, v_11=0.8661215489799035, v_12=-0.10948737831060629, v_22=0.18658686238607547),
t3=SymTensor3d_3(v_000=1.4639056597684772, v_001=0.39585204065168006, v_002=4.579797622976988, v_011=-2.3717172864430887, v_012=0.7010587169345798, v_022=0.9078116266746078, v_111=-0.6694856124915215, v_112=-1.4292770936051242, v_122=0.27363357183984083, v_222=-3.1505205293718657),
t4=SymTensor3d_4(v_0000=-9.4447716731816, v_0001=-4.973612689148469, v_0002=-20.92991278402753, v_0011=9.746922887744212, v_0012=-5.73747811922247, v_0022=-0.30215121456260485, v_0111=5.482856309197803, v_0112=5.262184073843096, v_0122=-0.5092436200493328, v_0222=15.667728710184438, v_1111=-8.358385283743033, v_1112=3.106386082300924, v_1122=-1.388537604001178, v_1222=2.631092036921547, v_2222=1.6906888185637832),
t5=SymTensor3d_5(v_00000=48.09523078243683, v_00001=41.53164634624548, v_00002=41.53164634624548, v_00011=-24.04761539121845, v_00012=15.843134845979208, v_00022=-24.04761539121845, v_00111=-33.327165801006274, v_00112=-8.204480545239235, v_00122=-8.204480545239235, v_00222=-33.327165801006274, v_01111=17.54061357947698, v_01112=-7.921567422989607, v_01122=6.507001811741461, v_01222=-7.921567422989607, v_02222=17.54061357947698, v_11111=28.630807971662435, v_11112=3.508122715895396, v_11122=4.696357829343838, v_11222=4.696357829343838, v_12222=3.508122715895396, v_22222=28.630807971662435)
)
Weirdly we can see that for dMk are different even though they will be contracted with the same a_k This is normal because we translate the moment dMk into the box A’, so even if we estimate the force in A’ after the translation we will still get the same force as the one we had in A before the translation. Which is arguably what we want XD.
1380 delta = dM_k_offset - dMp_k
1381
1382 print("delta =", delta)
1383 print("sqdist_t1 =", np.sqrt(delta.t1.inner(delta.t1)))
1384 print("sqdist_t2 =", np.sqrt(delta.t2.inner(delta.t2)) / 2)
1385 print("sqdist_t3 =", np.sqrt(delta.t3.inner(delta.t3)) / 6)
1386 print("sqdist_t4 =", np.sqrt(delta.t4.inner(delta.t4)) / 24)
1387 print("sqdist_t5 =", np.sqrt(delta.t5.inner(delta.t5)) / 120)
1388 print("(norm) =", tensor_collect_norm(delta))
delta = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=0.005066905810600875, v_1=0.005030136862605464, v_2=-0.009062158019562117),
t2=SymTensor3d_2(v_00=0.09792661753274357, v_01=-0.015603922553730554, v_02=0.11424626330652621, v_11=-0.05052894881105319, v_12=0.03500888091973883, v_22=-0.04739766872169093),
t3=SymTensor3d_3(v_000=-1.321115845648969, v_001=-0.4774812576091832, v_002=0.09422196374907532, v_011=0.24709178008571486, v_012=-0.2952641895649878, v_022=1.074024065563251, v_111=0.37944302590002776, v_112=-0.20214993284787508, v_122=0.09803823170915502, v_222=0.10792796909879865),
t4=SymTensor3d_4(v_0000=0.720369124576786, v_0001=2.307127570495746, v_0002=-7.551084521597305, v_0011=2.7353630804053086, v_0012=-0.6814668135150397, v_0022=-3.455732204982089, v_0111=-1.2052051534134964, v_0112=2.5256245155215677, v_0122=-1.1019224170822468, v_0222=5.025460006075747, v_1111=-2.889949819631486, v_1112=0.5155428243065439, v_1122=0.15458673922617927, v_1222=0.16592398920849716, v_2222=3.301145465755911),
t5=SymTensor3d_5(v_00000=29.37462722486708, v_00001=14.964485151102192, v_00002=23.82020554948329, v_00011=-18.653766853729984, v_00012=6.323647066867808, v_00022=-10.72086037113715, v_00111=-12.950667928217126, v_00112=-6.082380447810651, v_00122=-2.0138172228851037, v_00222=-17.737825101672673, v_01111=14.796015990427707, v_01112=-3.473597150653714, v_01122=3.8577508633022672, v_01222=-2.850049916214103, v_02222=6.8631095078348725, v_11111=11.532926856234116, v_11112=3.009104972431892, v_11122=1.4177410719830075, v_11222=3.073275475378759, v_12222=0.5960761509020944, v_22222=14.664549626293912)
)
sqdist_t1 = 0.011536833158261293
sqdist_t2 = 0.10420168152823184
sqdist_t3 = 0.4387428684956316
sqdist_t4 = 1.0125236327407858
sqdist_t5 = 1.1484704299114894
(norm) = 2.7154754458344
1392 a_i = (x_i[0] - s_A[0], x_i[1] - s_A[1], x_i[2] - s_A[2])
1393 a_ip = (x_i[0] - s_Ap[0], x_i[1] - s_Ap[1], x_i[2] - s_Ap[2])
1394
1395 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_i)
1396 a_kp = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_ip)
1397
1398 print("a_k =", a_k)
1399 print("a_kp =", a_kp)
a_k = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=0.19999999999999996, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.03999999999999998, v_01=0.039999999999999994, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0.007999999999999995, v_001=0.007999999999999997, v_002=0, v_011=0.008, v_012=0, v_022=0, v_111=0.008000000000000002, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=0.0015999999999999986, v_0001=0.001599999999999999, v_0002=0, v_0011=0.0015999999999999996, v_0012=0, v_0022=0, v_0111=0.0016, v_0112=0, v_0122=0, v_0222=0, v_1111=0.0016000000000000005, v_1112=0, v_1122=0, v_1222=0, v_2222=0)
)
a_kp = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=0.09999999999999987, v_1=0.1, v_2=0),
t2=SymTensor3d_2(v_00=0.009999999999999974, v_01=0.009999999999999988, v_02=0, v_11=0.010000000000000002, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0.0009999999999999961, v_001=0.0009999999999999974, v_002=0, v_011=0.000999999999999999, v_012=0, v_022=0, v_111=0.0010000000000000002, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=9.999999999999948e-05, v_0001=9.999999999999961e-05, v_0002=0, v_0011=9.999999999999976e-05, v_0012=0, v_0022=0, v_0111=9.99999999999999e-05, v_0112=0, v_0122=0, v_0222=0, v_1111=0.00010000000000000003, v_1112=0, v_1122=0, v_1222=0, v_2222=0)
)
1402 result = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
1403 resultp = shamrock.phys.contract_grav_moment_to_force_5(a_kp, dMp_k)
1404 result_offset = shamrock.phys.contract_grav_moment_to_force_5(a_kp, dM_k_offset)
1405
1406 print("force_i =", -Gconst * np.array(result))
1407 print("force_ip =", -Gconst * np.array(resultp))
1408 print("force_ip_offset =", -Gconst * np.array(result_offset), "force_i translated to A'")
force_i = [-0.68591296 -0.13718259 -0.34118049]
force_ip = [-0.68056702 -0.13639473 -0.3390527 ]
force_ip_offset = [-0.68591296 -0.13718259 -0.34118049] force_i translated to A'
As expected the delta is almost null
1412 delta_f = np.linalg.norm(np.array(result_offset) - np.array(result))
1413 delta_f /= np.linalg.norm(np.array(result))
1414 print("delta_f =", delta_f)
delta_f = 3.566330011800396e-17
Let’s modify FMM in a box to add the translation of the box A
1419 def test_grav_moment_offset(x_i, x_j, s_A, s_Ap, s_B, m_j, order, do_print):
1420
1421 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1422 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
1423 a_i = (x_i[0] - s_A[0], x_i[1] - s_A[1], x_i[2] - s_A[2])
1424 a_ip = (x_i[0] - s_Ap[0], x_i[1] - s_Ap[1], x_i[2] - s_Ap[2])
1425
1426 if do_print:
1427 print("x_i =", x_i)
1428 print("x_j =", x_j)
1429 print("s_A =", s_A)
1430 print("s_Ap =", s_Ap)
1431 print("s_B =", s_B)
1432 print("b_j =", b_j)
1433 print("r =", r)
1434 print("a_i =", a_i)
1435
1436 # compute the tensor product of the displacment
1437 if order == 1:
1438 Q_n_B = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_j)
1439 elif order == 2:
1440 Q_n_B = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_j)
1441 elif order == 3:
1442 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
1443 elif order == 4:
1444 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
1445 elif order == 5:
1446 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
1447 else:
1448 raise ValueError("Invalid order")
1449
1450 # multiply by mass to get the mass moment
1451 Q_n_B *= m_j
1452
1453 if do_print:
1454 print("Q_n_B =", Q_n_B)
1455
1456 # green function gradients
1457 if order == 1:
1458 D_n = shamrock.phys.green_func_grav_cartesian_1_1(r)
1459 elif order == 2:
1460 D_n = shamrock.phys.green_func_grav_cartesian_1_2(r)
1461 elif order == 3:
1462 D_n = shamrock.phys.green_func_grav_cartesian_1_3(r)
1463 elif order == 4:
1464 D_n = shamrock.phys.green_func_grav_cartesian_1_4(r)
1465 elif order == 5:
1466 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1467 else:
1468 raise ValueError("Invalid order")
1469
1470 if do_print:
1471 print("D_n =", D_n)
1472
1473 if order == 1:
1474 dM_k = shamrock.phys.get_dM_mat_1(D_n, Q_n_B)
1475 elif order == 2:
1476 dM_k = shamrock.phys.get_dM_mat_2(D_n, Q_n_B)
1477 elif order == 3:
1478 dM_k = shamrock.phys.get_dM_mat_3(D_n, Q_n_B)
1479 elif order == 4:
1480 dM_k = shamrock.phys.get_dM_mat_4(D_n, Q_n_B)
1481 elif order == 5:
1482 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
1483 else:
1484 raise ValueError("Invalid order")
1485
1486 if do_print:
1487 print("dM_k =", dM_k)
1488
1489 if order == 5:
1490 dM_k_offset = shamrock.phys.offset_dM_mat_5(dM_k, s_A, s_Ap)
1491 elif order == 4:
1492 dM_k_offset = shamrock.phys.offset_dM_mat_4(dM_k, s_A, s_Ap)
1493 elif order == 3:
1494 dM_k_offset = shamrock.phys.offset_dM_mat_3(dM_k, s_A, s_Ap)
1495 elif order == 2:
1496 dM_k_offset = shamrock.phys.offset_dM_mat_2(dM_k, s_A, s_Ap)
1497 else:
1498 raise ValueError("Invalid order")
1499
1500 if do_print:
1501 print("dM_k_offset =", dM_k_offset)
1502
1503 if order == 1:
1504 a_k = shamrock.math.SymTensorCollection_f64_0_0.from_vec(a_i)
1505 elif order == 2:
1506 a_k = shamrock.math.SymTensorCollection_f64_0_1.from_vec(a_i)
1507 elif order == 3:
1508 a_k = shamrock.math.SymTensorCollection_f64_0_2.from_vec(a_i)
1509 elif order == 4:
1510 a_k = shamrock.math.SymTensorCollection_f64_0_3.from_vec(a_i)
1511 elif order == 5:
1512 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_i)
1513 else:
1514 raise ValueError("Invalid order")
1515
1516 if do_print:
1517 print("a_k =", a_k)
1518
1519 if order == 1:
1520 a_kp = shamrock.math.SymTensorCollection_f64_0_0.from_vec(a_ip)
1521 elif order == 2:
1522 a_kp = shamrock.math.SymTensorCollection_f64_0_1.from_vec(a_ip)
1523 elif order == 3:
1524 a_kp = shamrock.math.SymTensorCollection_f64_0_2.from_vec(a_ip)
1525 elif order == 4:
1526 a_kp = shamrock.math.SymTensorCollection_f64_0_3.from_vec(a_ip)
1527 elif order == 5:
1528 a_kp = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_ip)
1529 else:
1530 raise ValueError("Invalid order")
1531
1532 if do_print:
1533 print("a_kp =", a_kp)
1534
1535 if order == 1:
1536 result = shamrock.phys.contract_grav_moment_to_force_1(a_k, dM_k)
1537 elif order == 2:
1538 result = shamrock.phys.contract_grav_moment_to_force_2(a_k, dM_k)
1539 elif order == 3:
1540 result = shamrock.phys.contract_grav_moment_to_force_3(a_k, dM_k)
1541 elif order == 4:
1542 result = shamrock.phys.contract_grav_moment_to_force_4(a_k, dM_k)
1543 elif order == 5:
1544 result = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
1545 else:
1546 raise ValueError("Invalid order")
1547
1548 Gconst = 1 # let's just set the grav constant to 1
1549 force_i = -Gconst * np.array(result)
1550 if do_print:
1551 print("force_i =", force_i)
1552
1553 if order == 1:
1554 result_offset = shamrock.phys.contract_grav_moment_to_force_1(a_kp, dM_k_offset)
1555 elif order == 2:
1556 result_offset = shamrock.phys.contract_grav_moment_to_force_2(a_kp, dM_k_offset)
1557 elif order == 3:
1558 result_offset = shamrock.phys.contract_grav_moment_to_force_3(a_kp, dM_k_offset)
1559 elif order == 4:
1560 result_offset = shamrock.phys.contract_grav_moment_to_force_4(a_kp, dM_k_offset)
1561 elif order == 5:
1562 result_offset = shamrock.phys.contract_grav_moment_to_force_5(a_kp, dM_k_offset)
1563 else:
1564 raise ValueError("Invalid order")
1565
1566 force_i_offset = -Gconst * np.array(result_offset)
1567 if do_print:
1568 print("force_i_offset =", force_i_offset)
1569
1570 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
1571 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1572 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
1573
1574 angle = (b_A_size + b_B_size) / b_dist
1575
1576 delta_A = np.linalg.norm(np.array(s_A) - np.array(s_Ap))
1577
1578 if do_print:
1579 print("b_A_size =", b_A_size)
1580 print("b_B_size =", b_B_size)
1581 print("b_dist =", b_dist)
1582 print("angle =", angle)
1583
1584 return force_i, force_i_offset, angle, delta_A
Let test for many different parameters. For clarification a perfect result here is that the translated dMk contracted with the new displacment ak_p give the same result as the original expansion (which it does ;) ).
1591 plt.figure()
1592 for order in range(2, 6):
1593 print("--------------------------------")
1594 print(f"Running FMM order = {order}")
1595 print("--------------------------------")
1596
1597 # set seed
1598 rng = np.random.default_rng(111)
1599
1600 N = 50000
1601
1602 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
1603 x_i_all = []
1604 for i in range(N):
1605 x_i_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1606
1607 # same for x_j
1608 x_j_all = []
1609 for i in range(N):
1610 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1611
1612 box_scale_fact_all = np.linspace(0, 0.1, N).tolist()
1613
1614 # same for box_1_center
1615 s_A_all = []
1616 s_Ap_all = []
1617 for p, box_scale_fact in zip(x_i_all, box_scale_fact_all):
1618 s_A_all.append(
1619 (
1620 p[0] + box_scale_fact * rng.uniform(-1, 1),
1621 p[1] + box_scale_fact * rng.uniform(-1, 1),
1622 p[2] + box_scale_fact * rng.uniform(-1, 1),
1623 )
1624 )
1625 s_Ap_all.append(
1626 (
1627 p[0] + box_scale_fact * rng.uniform(-1, 1),
1628 p[1] + box_scale_fact * rng.uniform(-1, 1),
1629 p[2] + box_scale_fact * rng.uniform(-1, 1),
1630 )
1631 )
1632
1633 # same for box_2_center
1634 s_B_all = []
1635 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1636 s_B_all.append(
1637 (
1638 p[0] + box_scale_fact * rng.uniform(-1, 1),
1639 p[1] + box_scale_fact * rng.uniform(-1, 1),
1640 p[2] + box_scale_fact * rng.uniform(-1, 1),
1641 )
1642 )
1643
1644 angles = []
1645 delta_A_all = []
1646 rel_errors = []
1647
1648 for x_i, x_j, s_A, s_Ap, s_B in zip(x_i_all, x_j_all, s_A_all, s_Ap_all, s_B_all):
1649
1650 force_i, force_i_offset, angle, delta_A = test_grav_moment_offset(
1651 x_i, x_j, s_A, s_Ap, s_B, m_j, order, do_print=False
1652 )
1653
1654 abs_error = np.linalg.norm(force_i_offset - force_i)
1655 rel_error = abs_error / np.linalg.norm(force_i)
1656
1657 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
1658 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1659 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
1660 angle = (b_A_size + b_B_size) / b_dist
1661
1662 if angle > 5.0 or angle < 1e-4:
1663 continue
1664
1665 angles.append(angle)
1666 delta_A_all.append(delta_A)
1667 rel_errors.append(rel_error)
1668
1669 print(f"Computed for {len(angles)} cases")
1670
1671 plt.scatter(angles, rel_errors, s=1, label=f"FMM order = {order}")
1672
1673
1674 plt.xlabel("Angle")
1675 plt.ylabel("$|f_{\\rm fmm} - f_{\\rm fmm, offset}| / |f_{\\rm fmm}|$ (Relative error)")
1676 plt.xscale("log")
1677 plt.yscale("log")
1678 plt.title("Grav moment translation error")
1679 plt.legend(loc="lower right")
1680 plt.grid()
1681 plt.show()
--------------------------------
Running FMM order = 2
--------------------------------
Computed for 49965 cases
--------------------------------
Running FMM order = 3
--------------------------------
Computed for 49965 cases
--------------------------------
Running FMM order = 4
--------------------------------
Computed for 49965 cases
--------------------------------
Running FMM order = 5
--------------------------------
Computed for 49965 cases
Small note on multipole method (Without FMM)#
In some ways a MM method is a FMM where there is no box A. If we reuse the formula at the start of this page we get:
\[\mathbf{b}_j = \mathbf{x}_j - \mathbf{s}_b\]
and the distance between the centers of the boxes is:
\[\mathbf{r} = \mathbf{s}_b - \mathbf{x}_i\]
This implies that the distance between the two particles is:
\[\mathbf{x}_j - \mathbf{x}_i = \mathbf{r} + \mathbf{b}_j\]
\[\begin{split}\phi(\mathbf{x}_i) &= - \mathcal{G}\iiint_V \rho(\mathbf{x}_j) G(\mathbf{x}_j - \mathbf{x}_i) ~{\rm d}^3\mathbf{x}_j\\
&\simeq - \mathcal{G}\iiint_V \rho(\mathbf{x}_j) \sum_{n = 0}^p \frac{1}{n!} \nabla_r^{(n)} G(\mathbf{r}) \cdot \mathbf{b}_j^{(n)} ~{\rm d}^3\mathbf{x}_j\\
&= - \mathcal{G}\sum_{n = 0}^p \frac{1}{n!} \underbrace{\nabla_r^{(n)} G(\mathbf{r})}_{D_n} \cdot \underbrace{\left(\iiint_V \rho(\mathbf{x}_j) \mathbf{b}_j^{(n)}~{\rm d}^3\mathbf{x}_j\right)}_{Q_n^B},\end{split}\]
where \(D_n\) are the gradients of the Green function and \(Q^B_n\) are the moments of the mass distribution. Hence, the force can then be written as follows:
\[\begin{split}f_{\rm g}(\mathbf{x}_i) &= -\nabla_i \phi(\mathbf{x}_i)\\
&= \mathcal{G}\iiint_V \rho(\mathbf{x}_j) \nabla_i G(\mathbf{x}_j - \mathbf{x}_i) ~{\rm d}^3\mathbf{x}_j\\
&= - \mathcal{G}\iiint_V \rho(\mathbf{x}_j) \nabla_j G(\mathbf{x}_j - \mathbf{x}_i) ~{\rm d}^3\mathbf{x}_j\\
&= - \mathcal{G}\sum_{n = 0}^p \frac{1}{n!} \nabla_r {D_n} \cdot {Q_n^B}\\
&= -\mathcal{G} \sum_{n = 0}^p \frac{1}{n!} {D_{n+1}} \cdot {Q_n^B} \\
&= -\mathcal{G} \sum_{n = 0}^p \frac{1}{n!} {Q_n^B} \cdot {D_{n+1}}\end{split}\]
As we can see, the expression of the MM force is litteraly the same contraction as the end of the FMM. Essentially in MM the green function moments are the equivalent of \(dM_k\) in FMM. So we can use the same final function but put the mass moments instead of box a displacements.
Q_n_B = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=0.2, v_1=0.2, v_2=0),
t2=SymTensor3d_2(v_00=0.04000000000000001, v_01=0.04000000000000001, v_02=0, v_11=0.04000000000000001, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0.008000000000000002, v_001=0.008000000000000002, v_002=0, v_011=0.008000000000000002, v_012=0, v_022=0, v_111=0.008000000000000002, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=0.0016000000000000005, v_0001=0.0016000000000000005, v_0002=0, v_0011=0.0016000000000000005, v_0012=0, v_0022=0, v_0111=0.0016000000000000005, v_0112=0, v_0122=0, v_0222=0, v_1111=0.0016000000000000005, v_1112=0, v_1122=0, v_1222=0, v_2222=0)
)
1747 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1748 print("D_n =", D_n)
D_n = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=0.6664823809676648, v_1=0.11108039682794413, v_2=-0),
t2=SymTensor3d_2(v_00=1.0657713749708153, v_01=0.27019555985175603, v_02=-0, v_11=-0.5103693908310947, v_12=-0, v_22=-0.5554019841397206),
t3=SymTensor3d_3(v_000=2.5193910310501564, v_001=0.8702244382612861, v_002=0, v_011=-1.1684132317913773, v_012=-0, v_022=-1.3509777992587801, v_111=-0.6450614717181563, v_112=0, v_122=-0.22516296654313003, v_222=0),
t4=SymTensor3d_4(v_0000=7.818204247354045, v_0001=3.4785951368788894, v_0002=0, v_0011=-3.467082056047611, v_0012=-0, v_0022=-4.35112219130643, v_0111=-2.565772299541876, v_0112=0, v_0122=-0.9128228373370134, v_0222=0, v_1111=2.4934043628881306, v_1112=0, v_1122=0.9736776931594812, v_1222=0, v_2222=3.37744449814695),
t5=SymTensor3d_5(v_00000=29.815100928798046, v_00001=16.56450061106264, v_00002=-0, v_00011=-12.422125244403624, v_00012=-0, v_00022=-17.392975684394447, v_00111=-12.144299934768547, v_00112=0, v_00122=-4.420200676294097, v_00222=0, v_01111=8.72149212006438, v_01112=0, v_01122=3.7006331243392436, v_01222=0, v_02222=13.692342560055202, v_11111=10.006156351816982, v_11112=-0, v_11122=2.138143582951563, v_11222=0, v_12222=2.2820570933425337, v_22222=0)
)
force_i = [-1.00133257e+00 -5.70430926e-04 -0.00000000e+00]
We can check that this is equivalent to the FMM with s_A = (0,0,0)
1758 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
1759 print("dM_k =", dM_k)
1760 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec((0, 0, 0))
1761 print("a_k =", a_k)
1762 force_i_fmm_sA_null = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
1763 Gconst = 1 # let's just set the grav constant to 1
1764 force_i_fmm_sA_null = -Gconst * np.array(force_i_fmm_sA_null)
1765 print("force_i_fmm_sA_null =", force_i_fmm_sA_null)
dM_k = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=1.0013325726540552, v_1=0.000570430925742737, v_2=0),
t2=SymTensor3d_2(v_00=-2.009991287593064, v_01=-0.025579931908720897, v_02=0, v_11=1.012081300493465, v_12=0, v_22=0.9979099870995988),
t3=SymTensor3d_3(v_000=5.7891904460271375, v_001=0.46404605817727873, v_002=0, v_011=-2.9347687627868018, v_012=0, v_022=-2.8544216832403366, v_111=-0.3534382459053615, v_112=0, v_122=-0.11060781227191739, v_222=0),
t4=SymTensor3d_4(v_0000=-17.094124555326182, v_0001=-4.307070210210693, v_0002=0, v_0011=8.380367091882045, v_0012=0, v_0022=8.71375746344414, v_0111=3.2503338624827096, v_0112=-0, v_0122=1.056736347727984, v_0222=-0, v_1111=-6.238934057264403, v_1112=-0, v_1122=-2.1414330346176427, v_1222=-0, v_2222=-6.572324428826498),
t5=SymTensor3d_5(v_00000=29.815100928798046, v_00001=16.56450061106264, v_00002=-0, v_00011=-12.422125244403624, v_00012=-0, v_00022=-17.392975684394447, v_00111=-12.144299934768547, v_00112=0, v_00122=-4.420200676294097, v_00222=0, v_01111=8.72149212006438, v_01112=0, v_01122=3.7006331243392436, v_01222=0, v_02222=13.692342560055202, v_11111=10.006156351816982, v_11112=-0, v_11122=2.138143582951563, v_11222=0, v_12222=2.2820570933425337, v_22222=0)
)
a_k = SymTensorCollection_f64_0_4(
t0=1,
t1=SymTensor3d_1(v_0=0, v_1=0, v_2=0),
t2=SymTensor3d_2(v_00=0, v_01=0, v_02=0, v_11=0, v_12=0, v_22=0),
t3=SymTensor3d_3(v_000=0, v_001=0, v_002=0, v_011=0, v_012=0, v_022=0, v_111=0, v_112=0, v_122=0, v_222=0),
t4=SymTensor3d_4(v_0000=0, v_0001=0, v_0002=0, v_0011=0, v_0012=0, v_0022=0, v_0111=0, v_0112=0, v_0122=0, v_0222=0, v_1111=0, v_1112=0, v_1122=0, v_1222=0, v_2222=0)
)
force_i_fmm_sA_null = [-1.00133257e+00 -5.70430926e-04 -0.00000000e+00]
Now we just need the analytical force to compare
1770 def analytic_force_i(x_i, x_j, Gconst):
1771 force_i_direct = (x_j[0] - x_i[0], x_j[1] - x_i[1], x_j[2] - x_i[2])
1772 force_i_direct /= np.linalg.norm(force_i_direct) ** 3
1773 force_i_direct *= m_i
1774 return force_i_direct
1775
1776
1777 force_i_exact = analytic_force_i(x_i, x_j, Gconst)
1778 print("force_i_exact =", force_i_exact)
force_i_exact = [-1. 0. 0.]
1781 abs_error = np.linalg.norm(force_i - force_i_exact)
1782 rel_error = abs_error / np.linalg.norm(force_i)
1783
1784
1785 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1786 b_dist = np.linalg.norm(np.array(x_i) - np.array(s_B))
1787 angle = (b_B_size) / b_dist
1788
1789 print("abs_error =", abs_error)
1790 print("rel_error =", rel_error)
1791 print("angle =", angle)
1792
1793 assert rel_error < 1e-2
abs_error = 0.0014495314137262958
rel_error = 0.0014476021434907016
angle = 0.23249527748763862
Let’s code MM in a box
1799 def run_mm(x_i, x_j, s_B, m_j, order, do_print):
1800
1801 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1802 r = (s_B[0] - x_i[0], s_B[1] - x_i[1], s_B[2] - x_i[2])
1803
1804 if do_print:
1805 print("x_i =", x_i)
1806 print("x_j =", x_j)
1807 print("s_B =", s_B)
1808 print("b_j =", b_j)
1809 print("r =", r)
1810
1811 # compute the tensor product of the displacment
1812 if order == 1:
1813 Q_n_B = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_j)
1814 elif order == 2:
1815 Q_n_B = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_j)
1816 elif order == 3:
1817 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
1818 elif order == 4:
1819 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
1820 elif order == 5:
1821 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
1822 else:
1823 raise ValueError("Invalid order")
1824
1825 # multiply by mass to get the mass moment
1826 Q_n_B *= m_j
1827
1828 if do_print:
1829 print("Q_n_B =", Q_n_B)
1830
1831 # green function gradients
1832 if order == 1:
1833 D_n = shamrock.phys.green_func_grav_cartesian_1_1(r)
1834 elif order == 2:
1835 D_n = shamrock.phys.green_func_grav_cartesian_1_2(r)
1836 elif order == 3:
1837 D_n = shamrock.phys.green_func_grav_cartesian_1_3(r)
1838 elif order == 4:
1839 D_n = shamrock.phys.green_func_grav_cartesian_1_4(r)
1840 elif order == 5:
1841 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1842 else:
1843 raise ValueError("Invalid order")
1844
1845 if do_print:
1846 print("D_n =", D_n)
1847
1848 if order == 1:
1849 result = shamrock.phys.contract_grav_moment_to_force_1(Q_n_B, D_n)
1850 elif order == 2:
1851 result = shamrock.phys.contract_grav_moment_to_force_2(Q_n_B, D_n)
1852 elif order == 3:
1853 result = shamrock.phys.contract_grav_moment_to_force_3(Q_n_B, D_n)
1854 elif order == 4:
1855 result = shamrock.phys.contract_grav_moment_to_force_4(Q_n_B, D_n)
1856 elif order == 5:
1857 result = shamrock.phys.contract_grav_moment_to_force_5(Q_n_B, D_n)
1858 else:
1859 raise ValueError("Invalid order")
1860
1861 Gconst = 1 # let's just set the grav constant to 1
1862 force_i = -Gconst * np.array(result)
1863 if do_print:
1864 print("force_i =", force_i)
1865
1866 force_i_exact = analytic_force_i(x_i, x_j, Gconst)
1867 if do_print:
1868 print("force_i_exact =", force_i_exact)
1869
1870 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1871 b_dist = np.linalg.norm(np.array(x_i) - np.array(s_B))
1872
1873 angle = (b_B_size) / b_dist
1874
1875 if do_print:
1876 print("b_A_size =", b_A_size)
1877 print("b_B_size =", b_B_size)
1878 print("b_dist =", b_dist)
1879 print("angle =", angle)
1880
1881 return force_i, force_i_exact, angle
1886 plt.figure()
1887 for order in range(1, 6):
1888 print("--------------------------------")
1889 print(f"Running MM order = {order}")
1890 print("--------------------------------")
1891
1892 # set seed
1893 rng = np.random.default_rng(111)
1894
1895 N = 50000
1896
1897 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
1898 x_i_all = []
1899 for i in range(N):
1900 x_i_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1901
1902 # same for x_j
1903 x_j_all = []
1904 for i in range(N):
1905 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1906
1907 box_scale_fact_all = np.linspace(0, 0.1, N).tolist()
1908
1909 # same for box_2_center
1910 s_B_all = []
1911 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1912 s_B_all.append(
1913 (
1914 p[0] + box_scale_fact * rng.uniform(-1, 1),
1915 p[1] + box_scale_fact * rng.uniform(-1, 1),
1916 p[2] + box_scale_fact * rng.uniform(-1, 1),
1917 )
1918 )
1919
1920 angles = []
1921 rel_errors = []
1922
1923 for x_i, x_j, s_B in zip(x_i_all, x_j_all, s_B_all):
1924
1925 force_i, force_i_exact, angle = run_mm(x_i, x_j, s_B, m_j, order, do_print=False)
1926
1927 abs_error = np.linalg.norm(force_i - force_i_exact)
1928 rel_error = abs_error / np.linalg.norm(force_i_exact)
1929
1930 if angle > 5.0 or angle < 1e-4:
1931 continue
1932
1933 angles.append(angle)
1934 rel_errors.append(rel_error)
1935
1936 print(f"Computed for {len(angles)} cases")
1937
1938 plt.scatter(angles, rel_errors, s=1, label=f"MM order = {order}")
1939
1940
1941 def plot_powerlaw(order, center_y):
1942 X = [1e-3, 1e-2 / 3, 1e-1]
1943 Y = [center_y * (x) ** order for x in X]
1944 plt.plot(X, Y, linestyle="dashed", color="black")
1945 bbox = dict(boxstyle="round", fc="blanchedalmond", ec="orange", alpha=0.9)
1946 plt.text(X[1], Y[1], f"$\\propto x^{order}$", fontsize=9, bbox=bbox)
1947
1948
1949 plot_powerlaw(1, 1)
1950 plot_powerlaw(2, 1)
1951 plot_powerlaw(3, 1)
1952 plot_powerlaw(4, 1)
1953 plot_powerlaw(5, 1)
1954
1955 plt.xlabel("Angle")
1956 plt.ylabel("Relative Error")
1957 plt.xscale("log")
1958 plt.yscale("log")
1959 plt.title("MM precision")
1960 plt.legend(loc="lower right")
1961 plt.grid()
1962 plt.show()
--------------------------------
Running MM order = 1
--------------------------------
Computed for 49923 cases
--------------------------------
Running MM order = 2
--------------------------------
Computed for 49923 cases
--------------------------------
Running MM order = 3
--------------------------------
Computed for 49923 cases
--------------------------------
Running MM order = 4
--------------------------------
Computed for 49923 cases
--------------------------------
Running MM order = 5
--------------------------------
Computed for 49923 cases
Total running time of the script: (0 minutes 58.088 seconds)
Estimated memory usage: 147 MB