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.2600000000000002, 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.2600000000000002, 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 2.08166817e-17 -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 2.08166817e-17 -0.00000000e+00]
force_i_exact = [-1. 0. 0.]
abs error = 1.1295700899707568e-16
rel error = 1.1295700899707568e-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 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
578 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
579 a_i = (x_i[0] - s_A[0], x_i[1] - s_A[1], x_i[2] - s_A[2])
580
581 if do_print:
582 print("x_i =", x_i)
583 print("x_j =", x_j)
584 print("s_A =", s_A)
585 print("s_B =", s_B)
586 print("b_j =", b_j)
587 print("r =", r)
588 print("a_i =", a_i)
589
590 # compute the tensor product of the displacment
591 if order == 1:
592 Q_n_B = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_j)
593 elif order == 2:
594 Q_n_B = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_j)
595 elif order == 3:
596 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
597 elif order == 4:
598 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
599 elif order == 5:
600 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
601 else:
602 raise ValueError("Invalid order")
603
604 # multiply by mass to get the mass moment
605 Q_n_B *= m_j
606
607 if do_print:
608 print("Q_n_B =", Q_n_B)
609
610 # green function gradients
611 if order == 1:
612 D_n = shamrock.phys.green_func_grav_cartesian_1_1(r)
613 elif order == 2:
614 D_n = shamrock.phys.green_func_grav_cartesian_1_2(r)
615 elif order == 3:
616 D_n = shamrock.phys.green_func_grav_cartesian_1_3(r)
617 elif order == 4:
618 D_n = shamrock.phys.green_func_grav_cartesian_1_4(r)
619 elif order == 5:
620 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
621 else:
622 raise ValueError("Invalid order")
623
624 if do_print:
625 print("D_n =", D_n)
626
627 if order == 1:
628 dM_k = shamrock.phys.get_dM_mat_1(D_n, Q_n_B)
629 elif order == 2:
630 dM_k = shamrock.phys.get_dM_mat_2(D_n, Q_n_B)
631 elif order == 3:
632 dM_k = shamrock.phys.get_dM_mat_3(D_n, Q_n_B)
633 elif order == 4:
634 dM_k = shamrock.phys.get_dM_mat_4(D_n, Q_n_B)
635 elif order == 5:
636 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
637 else:
638 raise ValueError("Invalid order")
639
640 if do_print:
641 print("dM_k =", dM_k)
642
643 if order == 1:
644 a_k = shamrock.math.SymTensorCollection_f64_0_0.from_vec(a_i)
645 elif order == 2:
646 a_k = shamrock.math.SymTensorCollection_f64_0_1.from_vec(a_i)
647 elif order == 3:
648 a_k = shamrock.math.SymTensorCollection_f64_0_2.from_vec(a_i)
649 elif order == 4:
650 a_k = shamrock.math.SymTensorCollection_f64_0_3.from_vec(a_i)
651 elif order == 5:
652 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_i)
653 else:
654 raise ValueError("Invalid order")
655
656 if do_print:
657 print("a_k =", a_k)
658
659 if order == 1:
660 result = shamrock.phys.contract_grav_moment_to_force_1(a_k, dM_k)
661 elif order == 2:
662 result = shamrock.phys.contract_grav_moment_to_force_2(a_k, dM_k)
663 elif order == 3:
664 result = shamrock.phys.contract_grav_moment_to_force_3(a_k, dM_k)
665 elif order == 4:
666 result = shamrock.phys.contract_grav_moment_to_force_4(a_k, dM_k)
667 elif order == 5:
668 result = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
669 else:
670 raise ValueError("Invalid order")
671
672 Gconst = 1 # let's just set the grav constant to 1
673 force_i = -Gconst * np.array(result)
674 if do_print:
675 print("force_i =", force_i)
676
677 force_i_exact = analytic_force_i(x_i, x_j, Gconst)
678 if do_print:
679 print("force_i_exact =", force_i_exact)
680
681 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
682 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
683 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
684
685 angle = (b_A_size + b_B_size) / b_dist
686
687 if do_print:
688 print("b_A_size =", b_A_size)
689 print("b_B_size =", b_B_size)
690 print("b_dist =", b_dist)
691 print("angle =", angle)
692
693 return force_i, force_i_exact, angle
Let’s try with some new parameters
698 s_B = (0, 0, 0)
699 s_A = (1, 0, 0)
700
701 box_B_size = 0.5
702 box_A_size = 0.5
703
704 x_j = (0.2, 0.2, 0.0)
705 x_i = (1.2, 0.2, 0.2)
706 m_j = 1
707 m_i = 1
708
709 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order=5, do_print=True)
710 ax, rel_error, abs_error = plot_fmm_case(
711 s_A, box_A_size, x_i, s_B, box_B_size, x_j, force_i, force_i_exact, 0.5
712 )
713
714 plt.title(f"FMM angle={angle:.5f} rel error={rel_error:.2e}")
715 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
716 plt.show()
717
718 print("force_i =", force_i)
719 print("force_i_exact =", force_i_exact)
720 print("abs error =", abs_error)
721 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.2600000000000002, 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.21430643e-17 -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.21430643e-17 -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#
728 s_B = (0, 0, 0)
729 s_A = (1, 0, 0)
730
731 box_B_size = 0.5
732 box_A_size = 0.5
733
734 x_j = (0.2, 0.2, 0.0)
735 x_i = (0.8, 0.2, 0.2)
736 m_j = 1
737 m_i = 1
738
739
740 for order in range(1, 6):
741 print("--------------------------------")
742 print(f"Running FMM order = {order}")
743 print("--------------------------------")
744
745 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order, do_print=True)
746 ax, rel_error, abs_error = plot_fmm_case(
747 s_A, box_A_size, x_i, s_B, box_B_size, x_j, force_i, force_i_exact, 0.2
748 )
749
750 plt.title(f"FMM order={order} angle={angle:.5f} rel error={rel_error:.2e}")
751 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
752 plt.show()
753
754 print("force_i =", force_i)
755 print("force_i_exact =", force_i_exact)
756 print("abs error =", abs_error)
757 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 1.38777878e-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 1.38777878e-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 1.73472348e-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 1.73472348e-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.2600000000000002, 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 6.07153217e-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 6.07153217e-17 -7.24000000e-01]
force_i_exact = [-2.37170825 0. -0.79056942]
abs error = 0.07940820523782492
rel error = 0.03176328209512998
Sweeping through angles#
765 s_B = (0, 0, 0)
766 s_A_all = [(0.8, 0, 0), (1, 0, 0), (1.2, 0, 0)]
767
768 box_B_size = 0.5
769 box_A_size = 0.5
770
771 x_j = (0.2, 0.2, 0.0)
772 x_i = (0.8, 0.2, 0.2)
773 m_j = 1
774 m_i = 1
775
776 order = 3
777
778 for s_A in s_A_all:
779 print("--------------------------------")
780 print(f"Running FMM s_a = {s_A}")
781 print("--------------------------------")
782
783 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order, do_print=True)
784 ax, rel_error, abs_error = plot_fmm_case(
785 s_A, box_A_size, x_i, s_B, box_B_size, x_j, force_i, force_i_exact, 0.2
786 )
787
788 plt.title(f"FMM order={order} angle={angle:.5f} rel error={rel_error:.2e}")
789 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
790 plt.show()
791
792 print("force_i =", force_i)
793 print("force_i_exact =", force_i_exact)
794 print("abs error =", abs_error)
795 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.648437499999993, 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.648437499999993, 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 1.38777878e-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 1.38777878e-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.
809 plt.figure()
810 for order in range(1, 6):
811 print("--------------------------------")
812 print(f"Running FMM order = {order}")
813 print("--------------------------------")
814
815 # set seed
816 rng = np.random.default_rng(111)
817
818 N = 50000
819
820 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
821 x_i_all = []
822 for i in range(N):
823 x_i_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
824
825 # same for x_j
826 x_j_all = []
827 for i in range(N):
828 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
829
830 box_scale_fact_all = np.linspace(0, 0.1, N).tolist()
831
832 # same for box_1_center
833 s_A_all = []
834 for p, box_scale_fact in zip(x_i_all, box_scale_fact_all):
835 s_A_all.append(
836 (
837 p[0] + box_scale_fact * rng.uniform(-1, 1),
838 p[1] + box_scale_fact * rng.uniform(-1, 1),
839 p[2] + box_scale_fact * rng.uniform(-1, 1),
840 )
841 )
842
843 # same for box_2_center
844 s_B_all = []
845 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
846 s_B_all.append(
847 (
848 p[0] + box_scale_fact * rng.uniform(-1, 1),
849 p[1] + box_scale_fact * rng.uniform(-1, 1),
850 p[2] + box_scale_fact * rng.uniform(-1, 1),
851 )
852 )
853
854 angles = []
855 rel_errors = []
856
857 for x_i, x_j, s_A, s_B in zip(x_i_all, x_j_all, s_A_all, s_B_all):
858 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order, do_print=False)
859
860 abs_error = np.linalg.norm(force_i - force_i_exact)
861 rel_error = abs_error / np.linalg.norm(force_i_exact)
862
863 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
864 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
865 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
866 angle = (b_A_size + b_B_size) / b_dist
867
868 if angle > 5.0 or angle < 1e-4:
869 continue
870
871 angles.append(angle)
872 rel_errors.append(rel_error)
873
874 print(f"Computed for {len(angles)} cases")
875
876 plt.scatter(angles, rel_errors, s=1, label=f"FMM order = {order}")
877
878
879 def plot_powerlaw(order, center_y):
880 X = [1e-3, 1e-2 / 3, 1e-1]
881 Y = [center_y * (x) ** order for x in X]
882 plt.plot(X, Y, linestyle="dashed", color="black")
883 bbox = dict(boxstyle="round", fc="blanchedalmond", ec="orange", alpha=0.9)
884 plt.text(X[1], Y[1], f"$\\propto x^{order}$", fontsize=9, bbox=bbox)
885
886
887 plot_powerlaw(1, 1)
888 plot_powerlaw(2, 1)
889 plot_powerlaw(3, 1)
890 plot_powerlaw(4, 1)
891 plot_powerlaw(5, 1)
892
893 plt.xlabel("Angle")
894 plt.ylabel("Relative Error")
895 plt.xscale("log")
896 plt.yscale("log")
897 plt.title("FMM precision")
898 plt.legend(loc="lower right")
899 plt.grid()
900 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.
930 s_B = (0, 0, 0)
931 box_B_size = 1.0
932 x_j = (0.2, 0.2, 0.0)
933 m_j = 1
934
935 s_B_new = (0.3, 0.3, 0.3)
def plot_mass_moment_offset(s_B, s_B_new, box_B_size):
944 def plot_mass_moment_offset(s_B, s_B_new, box_B_size):
945 box_B = shamrock.math.AABB_f64_3(
946 (
947 s_B[0] - box_B_size / 2,
948 s_B[1] - box_B_size / 2,
949 s_B[2] - box_B_size / 2,
950 ),
951 (
952 s_B[0] + box_B_size / 2,
953 s_B[1] + box_B_size / 2,
954 s_B[2] + box_B_size / 2,
955 ),
956 )
957
958 box_B_new = shamrock.math.AABB_f64_3(
959 (
960 s_B_new[0] - box_B_size / 2,
961 s_B_new[1] - box_B_size / 2,
962 s_B_new[2] - box_B_size / 2,
963 ),
964 (
965 s_B_new[0] + box_B_size / 2,
966 s_B_new[1] + box_B_size / 2,
967 s_B_new[2] + box_B_size / 2,
968 ),
969 )
970
971 fig = plt.figure()
972 ax = fig.add_subplot(111, projection="3d")
973
974 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
975 draw_arrow(ax, s_B_new, x_j, "red", "$b_j' = x_j - s_B'$")
976
977 ax.scatter(s_B[0], s_B[1], s_B[2], color="black", label="s_B")
978 ax.scatter(s_B_new[0], s_B_new[1], s_B_new[2], color="red", label="s_B'")
979 ax.scatter(x_j[0], x_j[1], x_j[2], color="blue", label="$x_j$")
980
981 draw_aabb(ax, box_B, "blue", 0.2)
982 draw_aabb(ax, box_B_new, "red", 0.2)
983
984 center_view = (0.0, 0.0, 0.0)
985 view_size = 2.0
986 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
987 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
988 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
989 ax.set_xlabel("X")
990 ax.set_ylabel("Y")
991 ax.set_zlabel("Z")
992
993 return ax
1003 plot_mass_moment_offset(s_B, s_B_new, box_B_size)
1004
1005 plt.title("Mass moment offset illustration")
1006 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
1007 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’
1027 Q_n_B_offset = shamrock.phys.offset_multipole_5(Q_n_B, s_B, s_B_new)
1028 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.010000000000000002, v_01=0.010000000000000002, v_02=0.03, v_11=0.010000000000000002, v_12=0.03, v_22=0.09),
t3=SymTensor3d_3(v_000=-0.0010000000000000044, v_001=-0.0010000000000000044, v_002=-0.0030000000000000027, v_011=-0.0010000000000000044, v_012=-0.0030000000000000027, v_022=-0.009000000000000001, v_111=-0.0010000000000000044, v_112=-0.0030000000000000027, v_122=-0.009000000000000001, v_222=-0.027),
t4=SymTensor3d_4(v_0000=0.00010000000000000221, v_0001=0.00010000000000000221, v_0002=0.0003000000000000021, v_0011=0.00010000000000000221, v_0012=0.0003000000000000021, v_0022=0.0009000000000000013, v_0111=0.00010000000000000221, v_0112=0.0003000000000000021, v_0122=0.0009000000000000013, v_0222=0.0027, v_1111=0.00010000000000000221, v_1112=0.0003000000000000021, v_1122=0.0009000000000000013, v_1222=0.0027, v_2222=0.0081),
t5=SymTensor3d_5(v_00000=-1.0000000000001056e-05, v_00001=-1.0000000000001056e-05, v_00002=-3.000000000000111e-05, v_00011=-1.0000000000001056e-05, v_00012=-3.000000000000111e-05, v_00022=-9.000000000000078e-05, v_00111=-1.0000000000001056e-05, v_00112=-3.000000000000111e-05, v_00122=-9.000000000000078e-05, v_00222=-0.0002700000000000006, v_01111=-1.0000000000001056e-05, v_01112=-3.000000000000111e-05, v_01122=-9.000000000000078e-05, v_01222=-0.0002700000000000006, v_02222=-0.0008100000000000002, v_11111=-1.0000000000001056e-05, v_11112=-3.000000000000111e-05, v_11122=-9.000000000000078e-05, v_11222=-0.0002700000000000006, v_12222=-0.0008100000000000002, v_22222=-0.00243)
)
Print the norm of the moment in box B’
1034 def tensor_collect_norm(d):
1035 # detect the type of the tensor collection
1036 if isinstance(d, shamrock.math.SymTensorCollection_f64_0_5):
1037 return (
1038 np.sqrt(d.t0 * d.t0)
1039 + np.sqrt(d.t1.inner(d.t1))
1040 + np.sqrt(d.t2.inner(d.t2)) / 2
1041 + np.sqrt(d.t3.inner(d.t3)) / 6
1042 + np.sqrt(d.t4.inner(d.t4)) / 24
1043 + np.sqrt(d.t5.inner(d.t5)) / 120
1044 )
1045 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_4):
1046 return (
1047 np.sqrt(d.t0 * d.t0)
1048 + np.sqrt(d.t1.inner(d.t1))
1049 + np.sqrt(d.t2.inner(d.t2)) / 2
1050 + np.sqrt(d.t3.inner(d.t3)) / 6
1051 + np.sqrt(d.t4.inner(d.t4)) / 24
1052 )
1053 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_3):
1054 return (
1055 np.sqrt(d.t0 * d.t0)
1056 + np.sqrt(d.t1.inner(d.t1))
1057 + np.sqrt(d.t2.inner(d.t2)) / 2
1058 + np.sqrt(d.t3.inner(d.t3)) / 6
1059 )
1060 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_2):
1061 return np.sqrt(d.t0 * d.t0) + np.sqrt(d.t1.inner(d.t1)) + np.sqrt(d.t2.inner(d.t2)) / 2
1062 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_1):
1063 return np.sqrt(d.t0 * d.t0) + np.sqrt(d.t1.inner(d.t1))
1064 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_0):
1065 return np.sqrt(d.t0 * d.t0)
1066 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_5):
1067 return (
1068 np.sqrt(d.t1.inner(d.t1))
1069 + np.sqrt(d.t2.inner(d.t2)) / 2
1070 + np.sqrt(d.t3.inner(d.t3)) / 6
1071 + np.sqrt(d.t4.inner(d.t4)) / 24
1072 + np.sqrt(d.t5.inner(d.t5)) / 120
1073 )
1074 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_4):
1075 return (
1076 np.sqrt(d.t1.inner(d.t1))
1077 + np.sqrt(d.t2.inner(d.t2)) / 2
1078 + np.sqrt(d.t3.inner(d.t3)) / 6
1079 + np.sqrt(d.t4.inner(d.t4)) / 24
1080 )
1081 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_3):
1082 return (
1083 np.sqrt(d.t1.inner(d.t1))
1084 + np.sqrt(d.t2.inner(d.t2)) / 2
1085 + np.sqrt(d.t3.inner(d.t3)) / 6
1086 )
1087 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_2):
1088 return np.sqrt(d.t1.inner(d.t1)) + np.sqrt(d.t2.inner(d.t2)) / 2
1089 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_1):
1090 return np.sqrt(d.t1.inner(d.t1))
1091 else:
1092 raise ValueError(f"Unsupported tensor collection type: {type(d)}")
1093
1094
1095 print("Q_n_B norm =", tensor_collect_norm(Q_n_B))
1096 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
1100 delta = Q_n_B_offset - Q_n_Bp
1101 print("delta =", delta)
1102
1103
1104 sqdist_t0 = delta.t0 * delta.t0
1105 sqdist_t1 = delta.t1.inner(delta.t1)
1106 sqdist_t2 = delta.t2.inner(delta.t2)
1107 sqdist_t3 = delta.t3.inner(delta.t3)
1108 sqdist_t4 = delta.t4.inner(delta.t4)
1109 sqdist_t5 = delta.t5.inner(delta.t5)
1110 print("sqdist_t0 =", sqdist_t0)
1111 print("sqdist_t1 =", sqdist_t1)
1112 print("sqdist_t2 =", sqdist_t2)
1113 print("sqdist_t3 =", sqdist_t3)
1114 print("sqdist_t4 =", sqdist_t4)
1115 print("sqdist_t5 =", sqdist_t5)
1116
1117 norm_delta = (
1118 np.sqrt(sqdist_t0)
1119 + np.sqrt(sqdist_t1)
1120 + np.sqrt(sqdist_t2) / 2
1121 + np.sqrt(sqdist_t3) / 6
1122 + np.sqrt(sqdist_t4) / 24
1123 + np.sqrt(sqdist_t5) / 120
1124 )
1125 print("norm_delta =", norm_delta)
1126
1127 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=6.938893903907228e-18, v_01=6.938893903907228e-18, v_02=6.938893903907228e-18, v_11=6.938893903907228e-18, v_12=6.938893903907228e-18, v_22=0),
t3=SymTensor3d_3(v_000=-4.9873299934333204e-18, v_001=-4.9873299934333204e-18, v_002=-4.336808689942018e-18, v_011=-4.9873299934333204e-18, v_012=-4.336808689942018e-18, v_022=-3.469446951953614e-18, v_111=-4.9873299934333204e-18, v_112=-4.336808689942018e-18, v_122=-3.469446951953614e-18, v_222=0),
t4=SymTensor3d_4(v_0000=2.303929616531697e-18, v_0001=2.303929616531697e-18, v_0002=2.3310346708438345e-18, v_0011=2.303929616531697e-18, v_0012=2.3310346708438345e-18, v_0022=1.734723475976807e-18, v_0111=2.303929616531697e-18, v_0112=2.3310346708438345e-18, v_0122=1.734723475976807e-18, v_0222=8.673617379884035e-19, v_1111=2.303929616531697e-18, v_1112=2.3310346708438345e-18, v_1122=1.734723475976807e-18, v_1222=8.673617379884035e-19, v_2222=0),
t5=SymTensor3d_5(v_00000=-1.0672615135404184e-18, v_00001=-1.0672615135404184e-18, v_00002=-1.1384122811097797e-18, v_00011=-1.0672615135404184e-18, v_00012=-1.1384122811097797e-18, v_00022=-8.402566836762659e-19, v_00111=-1.0672615135404184e-18, v_00112=-1.1384122811097797e-18, v_00122=-8.402566836762659e-19, v_00222=-7.589415207398531e-19, v_01111=-1.0672615135404184e-18, v_01112=-1.1384122811097797e-18, v_01122=-8.402566836762659e-19, v_01222=-7.589415207398531e-19, v_02222=-4.336808689942018e-19, v_11111=-1.0672615135404184e-18, v_11112=-1.1384122811097797e-18, v_11122=-8.402566836762659e-19, v_11222=-7.589415207398531e-19, v_12222=-4.336808689942018e-19, v_22222=0)
)
sqdist_t0 = 0.0
sqdist_t1 = 0.0
sqdist_t2 = 3.851859888774472e-34
sqdist_t3 = 4.969049719795974e-34
sqdist_t4 = 3.370494952112745e-34
sqdist_t5 = 2.215310939620759e-34
norm_delta = 1.44172925989326e-17
rel error = 1.0347731059478237e-17
We now want to explore the precision of the offset as a function of the order & distance
1132 plt.figure()
1133
1134 for order in range(0, 6):
1135 # set seed
1136 rng = np.random.default_rng(111)
1137
1138 N = 50000
1139
1140 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
1141 x_j_all = []
1142 for i in range(N):
1143 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1144
1145 box_scale_fact_all = np.linspace(0, 1, N).tolist()
1146
1147 # same for box_1_center
1148 s_B_all = []
1149 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1150 s_B_all.append(
1151 (
1152 p[0] + box_scale_fact * rng.uniform(-1, 1),
1153 p[1] + box_scale_fact * rng.uniform(-1, 1),
1154 p[2] + box_scale_fact * rng.uniform(-1, 1),
1155 )
1156 )
1157
1158 # same for box_2_center
1159 s_Bp_all = []
1160 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1161 s_Bp_all.append(
1162 (
1163 p[0] + box_scale_fact * rng.uniform(-1, 1),
1164 p[1] + box_scale_fact * rng.uniform(-1, 1),
1165 p[2] + box_scale_fact * rng.uniform(-1, 1),
1166 )
1167 )
1168
1169 center_distances = []
1170 rel_errors = []
1171 for x_j, s_B, s_Bp in zip(x_j_all, s_B_all, s_Bp_all):
1172 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1173 b_jp = (x_j[0] - s_Bp[0], x_j[1] - s_Bp[1], x_j[2] - s_Bp[2])
1174
1175 if order == 5:
1176 Q_n_B = shamrock.math.SymTensorCollection_f64_0_5.from_vec(b_j)
1177 Q_n_B *= m_j
1178
1179 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_5.from_vec(b_jp)
1180 Q_n_Bp *= m_j
1181
1182 Q_n_B_offset = shamrock.phys.offset_multipole_5(Q_n_B, s_B, s_Bp)
1183 elif order == 4:
1184 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
1185 Q_n_B *= m_j
1186
1187 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_jp)
1188 Q_n_Bp *= m_j
1189
1190 Q_n_B_offset = shamrock.phys.offset_multipole_4(Q_n_B, s_B, s_Bp)
1191 elif order == 3:
1192 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
1193 Q_n_B *= m_j
1194
1195 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_jp)
1196 Q_n_Bp *= m_j
1197
1198 Q_n_B_offset = shamrock.phys.offset_multipole_3(Q_n_B, s_B, s_Bp)
1199 elif order == 2:
1200 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
1201 Q_n_B *= m_j
1202
1203 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_jp)
1204 Q_n_Bp *= m_j
1205
1206 Q_n_B_offset = shamrock.phys.offset_multipole_2(Q_n_B, s_B, s_Bp)
1207 elif order == 1:
1208 Q_n_B = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_j)
1209 Q_n_B *= m_j
1210
1211 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_jp)
1212 Q_n_Bp *= m_j
1213
1214 Q_n_B_offset = shamrock.phys.offset_multipole_1(Q_n_B, s_B, s_Bp)
1215 elif order == 0:
1216 Q_n_B = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_j)
1217 Q_n_B *= m_j
1218
1219 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_jp)
1220 Q_n_Bp *= m_j
1221
1222 Q_n_B_offset = shamrock.phys.offset_multipole_0(Q_n_B, s_B, s_Bp)
1223 else:
1224 raise ValueError(f"Unsupported offset order: {order}")
1225
1226 delta = Q_n_B_offset - Q_n_Bp
1227
1228 rel_error = tensor_collect_norm(delta) / tensor_collect_norm(Q_n_B)
1229 rel_errors.append(rel_error)
1230
1231 center_distances.append(np.linalg.norm(np.array(s_B) - np.array(s_Bp)))
1232
1233 plt.scatter(center_distances, rel_errors, s=1, label=f"multipole order = {order}")
1234
1235 plt.xlabel("$\\vert \\vert s_B - s_B'\\vert \\vert$")
1236 plt.ylabel(
1237 "$\\vert \\vert Q_n(s_B) - Q_n(s_B') \\vert \\vert / \\vert \\vert Q_n(s_B) \\vert \\vert$ (relative error) "
1238 )
1239 plt.xscale("log")
1240 plt.yscale("log")
1241 plt.title("Mass moment offset precision")
1242 plt.legend(loc="lower right")
1243 plt.grid()
1244 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.
1273 def plot_grav_moment_offset(s_A, s_Ap, s_B, box_A_size, box_B_size, x_j):
1274 box_A = shamrock.math.AABB_f64_3(
1275 (
1276 s_A[0] - box_A_size / 2,
1277 s_A[1] - box_A_size / 2,
1278 s_A[2] - box_A_size / 2,
1279 ),
1280 (
1281 s_A[0] + box_A_size / 2,
1282 s_A[1] + box_A_size / 2,
1283 s_A[2] + box_A_size / 2,
1284 ),
1285 )
1286
1287 box_Ap = shamrock.math.AABB_f64_3(
1288 (
1289 s_Ap[0] - box_A_size / 2,
1290 s_Ap[1] - box_A_size / 2,
1291 s_Ap[2] - box_A_size / 2,
1292 ),
1293 (
1294 s_Ap[0] + box_A_size / 2,
1295 s_Ap[1] + box_A_size / 2,
1296 s_Ap[2] + box_A_size / 2,
1297 ),
1298 )
1299
1300 box_B = shamrock.math.AABB_f64_3(
1301 (
1302 s_B[0] - box_B_size / 2,
1303 s_B[1] - box_B_size / 2,
1304 s_B[2] - box_B_size / 2,
1305 ),
1306 (
1307 s_B[0] + box_B_size / 2,
1308 s_B[1] + box_B_size / 2,
1309 s_B[2] + box_B_size / 2,
1310 ),
1311 )
1312
1313 fig = plt.figure()
1314 ax = fig.add_subplot(111, projection="3d")
1315
1316 draw_arrow(ax, s_A, s_B, "purple", "$r = s_B - s_A$")
1317 draw_arrow(ax, s_Ap, s_B, "purple", "$r' = s_B - s_A'$")
1318
1319 ax.scatter(s_A[0], s_A[1], s_A[2], color="black", label="s_A")
1320 ax.scatter(s_Ap[0], s_Ap[1], s_Ap[2], color="black", label="s_Ap")
1321 ax.scatter(s_B[0], s_B[1], s_B[2], color="black", label="s_B")
1322
1323 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
1324
1325 ax.scatter(x_j[0], x_j[1], x_j[2], color="red", label="$x_j$")
1326
1327 draw_aabb(ax, box_A, "blue", 0.1)
1328 draw_aabb(ax, box_Ap, "cyan", 0.1)
1329 draw_aabb(ax, box_B, "red", 0.1)
1330
1331 center_view = (0.5, 0.0, 0.0)
1332 view_size = 2.0
1333 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
1334 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
1335 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
1336 ax.set_xlabel("X")
1337 ax.set_ylabel("Y")
1338 ax.set_zlabel("Z")
1339
1340 return ax, rel_error, abs_error
1341
1342
1343 s_B = (0, -0.2, -0.2)
1344 s_A = (1, 0, 0)
1345 s_Ap = (1.1, 0.1, 0.0)
1346
1347 box_B_size = 0.5
1348 box_A_size = 0.5
1349
1350 x_j = (0.2, 0.0, -0.5)
1351 x_i = (1.2, 0.2, 0.0)
1352 m_j = 1
1353 m_i = 1
1354
1355 plot_grav_moment_offset(s_A, s_Ap, s_B, box_A_size, box_B_size, x_j)
1356
1357 plt.title("Mass moment offset illustration")
1358 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
1359 plt.show()
1360
1361 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1362 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
1363 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)
)
1372 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1373 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
1374 # print("D_n =",D_n)
1375 print("dM_k =", dM_k)
dM_k = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=0.9330692614896893, v_1=-0.0067602090561549, v_2=0.597533732675385),
t2=SymTensor3d_2(v_00=-1.298112911867752, v_01=0.11610754537124779, v_02=-1.813416530995671, v_11=1.2391425006660395, v_12=0.007638654300740011, v_22=0.05897041120171276),
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.2300083683889491, v_112=-2.368831572596142, v_122=0.5109410987828298, v_222=-5.202772318170674),
t4=SymTensor3d_4(v_0000=-18.407459386049837, v_0001=-6.722015784651175, v_0002=-26.667390903250002, v_0011=15.484401006966687, v_0012=-6.501343549296468, v_0022=2.923058379083166, v_0111=7.061511531350737, v_0112=6.87478887066598, v_0122=-0.33949574669955473, v_0222=19.792602032584032, v_1111=-12.975527438856975, v_1112=3.5477305530103456, v_1122=-2.508873568109709, v_1222=2.9536129962861244, v_2222=-0.41418481097345694),
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.507001811741462, v_01222=-7.921567422989607, v_02222=17.54061357947698, v_11111=28.63080797166243, v_11112=3.508122715895396, v_11122=4.696357829343838, v_11222=4.696357829343838, v_12222=3.508122715895396, v_22222=28.63080797166243)
)
1378 Dp_n = shamrock.phys.green_func_grav_cartesian_1_5(rp)
1379 dMp_k = shamrock.phys.get_dM_mat_5(Dp_n, Q_n_B)
1380 # print("Dp_n =",Dp_n)
1381 print("dMp_k =", dMp_k)
dMp_k = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=0.8051957104860917, v_1=0.08591568124200472, v_2=0.45554150189781367),
t2=SymTensor3d_2(v_00=-1.1506350288987228, v_01=-0.18473251757367654, v_02=-1.2544427285969801, v_11=0.9166504977909566, v_12=-0.1444962592303451, v_22=0.2339845311077664),
t3=SymTensor3d_3(v_000=2.785021505417446, v_001=0.873333298260863, v_002=4.485575659227913, v_011=-2.618809066528804, v_012=0.9963229064995676, v_022=-0.166212438888643, v_111=-1.0489286383915493, v_112=-1.2271271607572491, v_122=0.17559534013068578, v_222=-3.2584484984706643),
t4=SymTensor3d_4(v_0000=-10.16514079775839, v_0001=-7.280740259644215, v_0002=-13.378828262430222, v_0011=7.011559807338903, v_0012=-5.056011305707431, v_0022=3.153580990419484, v_0111=6.688061462611301, v_0112=2.736559558321528, v_0122=0.5926787970329143, v_0222=10.642268704108691, v_1111=-5.468435464111547, v_1112=2.59084325799438, v_1122=-1.5431243432273574, v_1222=2.4651680477130493, v_2222=-1.6104566471921289),
t5=SymTensor3d_5(v_00000=18.720603557569746, v_00001=26.56716119514328, v_00002=17.711440796762187, v_00011=-5.39384853748847, 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.7445975890492718, v_01112=-4.447970272335893, v_01122=2.649250948439194, v_01222=-5.071517506775503, v_02222=10.677504071642106, v_11111=17.09788111542832, v_11112=0.4990177434635039, v_11122=3.2786167573608305, v_11222=1.6230823539650792, v_12222=2.9120465649933016, v_22222=13.966258345368521)
)
Offset the grav moment to dMp_k
1385 dM_k_offset = shamrock.phys.offset_dM_mat_5(dM_k, s_A, s_Ap)
1386 print("dM_k_offset =", dM_k_offset)
dM_k_offset = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=0.8102626162966925, v_1=0.09094581810461014, v_2=0.4464793438782518),
t2=SymTensor3d_2(v_00=-1.0527084113659793, v_01=-0.20033644012740717, v_02=-1.1401964652904537, v_11=0.8661215489799033, v_12=-0.10948737831060629, v_22=0.18658686238607528),
t3=SymTensor3d_3(v_000=1.4639056597684772, v_001=0.39585204065168017, v_002=4.579797622976988, v_011=-2.3717172864430887, v_012=0.7010587169345797, v_022=0.9078116266746078, v_111=-0.6694856124915218, v_112=-1.4292770936051244, v_122=0.273633571839841, v_222=-3.150520529371865),
t4=SymTensor3d_4(v_0000=-9.444771673181602, v_0001=-4.973612689148469, v_0002=-20.92991278402753, v_0011=9.746922887744212, v_0012=-5.7374781192224695, v_0022=-0.30215121456260485, v_0111=5.482856309197805, v_0112=5.262184073843095, v_0122=-0.5092436200493327, v_0222=15.667728710184441, v_1111=-8.358385283743033, v_1112=3.106386082300924, v_1122=-1.3885376040011783, v_1222=2.6310920369215465, v_2222=1.6906888185637818),
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.507001811741462, v_01222=-7.921567422989607, v_02222=17.54061357947698, v_11111=28.63080797166243, v_11112=3.508122715895396, v_11122=4.696357829343838, v_11222=4.696357829343838, v_12222=3.508122715895396, v_22222=28.63080797166243)
)
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.
1393 delta = dM_k_offset - dMp_k
1394
1395 print("delta =", delta)
1396 print("sqdist_t1 =", np.sqrt(delta.t1.inner(delta.t1)))
1397 print("sqdist_t2 =", np.sqrt(delta.t2.inner(delta.t2)) / 2)
1398 print("sqdist_t3 =", np.sqrt(delta.t3.inner(delta.t3)) / 6)
1399 print("sqdist_t4 =", np.sqrt(delta.t4.inner(delta.t4)) / 24)
1400 print("sqdist_t5 =", np.sqrt(delta.t5.inner(delta.t5)) / 120)
1401 print("(norm) =", tensor_collect_norm(delta))
delta = SymTensorCollection_f64_1_5(
t1=SymTensor3d_1(v_0=0.005066905810600764, v_1=0.005030136862605422, v_2=-0.009062158019561894),
t2=SymTensor3d_2(v_00=0.09792661753274357, v_01=-0.015603922553730637, v_02=0.11424626330652643, v_11=-0.0505289488110533, v_12=0.035008880919738805, v_22=-0.04739766872169113),
t3=SymTensor3d_3(v_000=-1.321115845648969, v_001=-0.47748125760918286, v_002=0.09422196374907532, v_011=0.2470917800857153, v_012=-0.2952641895649879, v_022=1.074024065563251, v_111=0.37944302590002743, v_112=-0.2021499328478753, v_122=0.09803823170915521, v_222=0.10792796909879909),
t4=SymTensor3d_4(v_0000=0.7203691245767878, v_0001=2.307127570495746, v_0002=-7.551084521597307, v_0011=2.7353630804053086, v_0012=-0.6814668135150388, v_0022=-3.455732204982089, v_0111=-1.2052051534134955, v_0112=2.525624515521567, v_0122=-1.101922417082247, v_0222=5.02546000607575, v_1111=-2.889949819631486, v_1112=0.5155428243065439, v_1122=0.15458673922617905, v_1222=0.16592398920849716, v_2222=3.3011454657559107),
t5=SymTensor3d_5(v_00000=29.374627224867083, v_00001=14.964485151102199, v_00002=23.820205549483294, v_00011=-18.65376685372998, 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.857750863302268, v_01222=-2.850049916214104, v_02222=6.8631095078348725, v_11111=11.532926856234113, v_11112=3.009104972431892, v_11122=1.4177410719830075, v_11222=3.073275475378759, v_12222=0.5960761509020944, v_22222=14.66454962629391)
)
sqdist_t1 = 0.011536833158261052
sqdist_t2 = 0.104201681528232
sqdist_t3 = 0.4387428684956316
sqdist_t4 = 1.012523632740786
sqdist_t5 = 1.1484704299114894
(norm) = 2.7154754458344
1405 a_i = (x_i[0] - s_A[0], x_i[1] - s_A[1], x_i[2] - s_A[2])
1406 a_ip = (x_i[0] - s_Ap[0], x_i[1] - s_Ap[1], x_i[2] - s_Ap[2])
1407
1408 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_i)
1409 a_kp = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_ip)
1410
1411 print("a_k =", a_k)
1412 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)
)
1415 result = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
1416 resultp = shamrock.phys.contract_grav_moment_to_force_5(a_kp, dMp_k)
1417 result_offset = shamrock.phys.contract_grav_moment_to_force_5(a_kp, dM_k_offset)
1418
1419 print("force_i =", -Gconst * np.array(result))
1420 print("force_ip =", -Gconst * np.array(resultp))
1421 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
1425 delta_f = np.linalg.norm(np.array(result_offset) - np.array(result))
1426 delta_f /= np.linalg.norm(np.array(result))
1427 print("delta_f =", delta_f)
delta_f = 1.6342977232268068e-16
Let’s modify FMM in a box to add the translation of the box A
1432 def test_grav_moment_offset(x_i, x_j, s_A, s_Ap, s_B, m_j, order, do_print):
1433 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1434 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
1435 a_i = (x_i[0] - s_A[0], x_i[1] - s_A[1], x_i[2] - s_A[2])
1436 a_ip = (x_i[0] - s_Ap[0], x_i[1] - s_Ap[1], x_i[2] - s_Ap[2])
1437
1438 if do_print:
1439 print("x_i =", x_i)
1440 print("x_j =", x_j)
1441 print("s_A =", s_A)
1442 print("s_Ap =", s_Ap)
1443 print("s_B =", s_B)
1444 print("b_j =", b_j)
1445 print("r =", r)
1446 print("a_i =", a_i)
1447
1448 # compute the tensor product of the displacment
1449 if order == 1:
1450 Q_n_B = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_j)
1451 elif order == 2:
1452 Q_n_B = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_j)
1453 elif order == 3:
1454 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
1455 elif order == 4:
1456 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
1457 elif order == 5:
1458 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
1459 else:
1460 raise ValueError("Invalid order")
1461
1462 # multiply by mass to get the mass moment
1463 Q_n_B *= m_j
1464
1465 if do_print:
1466 print("Q_n_B =", Q_n_B)
1467
1468 # green function gradients
1469 if order == 1:
1470 D_n = shamrock.phys.green_func_grav_cartesian_1_1(r)
1471 elif order == 2:
1472 D_n = shamrock.phys.green_func_grav_cartesian_1_2(r)
1473 elif order == 3:
1474 D_n = shamrock.phys.green_func_grav_cartesian_1_3(r)
1475 elif order == 4:
1476 D_n = shamrock.phys.green_func_grav_cartesian_1_4(r)
1477 elif order == 5:
1478 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1479 else:
1480 raise ValueError("Invalid order")
1481
1482 if do_print:
1483 print("D_n =", D_n)
1484
1485 if order == 1:
1486 dM_k = shamrock.phys.get_dM_mat_1(D_n, Q_n_B)
1487 elif order == 2:
1488 dM_k = shamrock.phys.get_dM_mat_2(D_n, Q_n_B)
1489 elif order == 3:
1490 dM_k = shamrock.phys.get_dM_mat_3(D_n, Q_n_B)
1491 elif order == 4:
1492 dM_k = shamrock.phys.get_dM_mat_4(D_n, Q_n_B)
1493 elif order == 5:
1494 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
1495 else:
1496 raise ValueError("Invalid order")
1497
1498 if do_print:
1499 print("dM_k =", dM_k)
1500
1501 if order == 5:
1502 dM_k_offset = shamrock.phys.offset_dM_mat_5(dM_k, s_A, s_Ap)
1503 elif order == 4:
1504 dM_k_offset = shamrock.phys.offset_dM_mat_4(dM_k, s_A, s_Ap)
1505 elif order == 3:
1506 dM_k_offset = shamrock.phys.offset_dM_mat_3(dM_k, s_A, s_Ap)
1507 elif order == 2:
1508 dM_k_offset = shamrock.phys.offset_dM_mat_2(dM_k, s_A, s_Ap)
1509 elif order == 1:
1510 dM_k_offset = shamrock.phys.offset_dM_mat_1(dM_k, s_A, s_Ap)
1511 else:
1512 raise ValueError("Invalid order")
1513
1514 if do_print:
1515 print("dM_k_offset =", dM_k_offset)
1516
1517 if order == 1:
1518 a_k = shamrock.math.SymTensorCollection_f64_0_0.from_vec(a_i)
1519 elif order == 2:
1520 a_k = shamrock.math.SymTensorCollection_f64_0_1.from_vec(a_i)
1521 elif order == 3:
1522 a_k = shamrock.math.SymTensorCollection_f64_0_2.from_vec(a_i)
1523 elif order == 4:
1524 a_k = shamrock.math.SymTensorCollection_f64_0_3.from_vec(a_i)
1525 elif order == 5:
1526 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_i)
1527 else:
1528 raise ValueError("Invalid order")
1529
1530 if do_print:
1531 print("a_k =", a_k)
1532
1533 if order == 1:
1534 a_kp = shamrock.math.SymTensorCollection_f64_0_0.from_vec(a_ip)
1535 elif order == 2:
1536 a_kp = shamrock.math.SymTensorCollection_f64_0_1.from_vec(a_ip)
1537 elif order == 3:
1538 a_kp = shamrock.math.SymTensorCollection_f64_0_2.from_vec(a_ip)
1539 elif order == 4:
1540 a_kp = shamrock.math.SymTensorCollection_f64_0_3.from_vec(a_ip)
1541 elif order == 5:
1542 a_kp = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_ip)
1543 else:
1544 raise ValueError("Invalid order")
1545
1546 if do_print:
1547 print("a_kp =", a_kp)
1548
1549 if order == 1:
1550 result = shamrock.phys.contract_grav_moment_to_force_1(a_k, dM_k)
1551 elif order == 2:
1552 result = shamrock.phys.contract_grav_moment_to_force_2(a_k, dM_k)
1553 elif order == 3:
1554 result = shamrock.phys.contract_grav_moment_to_force_3(a_k, dM_k)
1555 elif order == 4:
1556 result = shamrock.phys.contract_grav_moment_to_force_4(a_k, dM_k)
1557 elif order == 5:
1558 result = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
1559 else:
1560 raise ValueError("Invalid order")
1561
1562 Gconst = 1 # let's just set the grav constant to 1
1563 force_i = -Gconst * np.array(result)
1564 if do_print:
1565 print("force_i =", force_i)
1566
1567 if order == 1:
1568 result_offset = shamrock.phys.contract_grav_moment_to_force_1(a_kp, dM_k_offset)
1569 elif order == 2:
1570 result_offset = shamrock.phys.contract_grav_moment_to_force_2(a_kp, dM_k_offset)
1571 elif order == 3:
1572 result_offset = shamrock.phys.contract_grav_moment_to_force_3(a_kp, dM_k_offset)
1573 elif order == 4:
1574 result_offset = shamrock.phys.contract_grav_moment_to_force_4(a_kp, dM_k_offset)
1575 elif order == 5:
1576 result_offset = shamrock.phys.contract_grav_moment_to_force_5(a_kp, dM_k_offset)
1577 else:
1578 raise ValueError("Invalid order")
1579
1580 force_i_offset = -Gconst * np.array(result_offset)
1581 if do_print:
1582 print("force_i_offset =", force_i_offset)
1583
1584 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
1585 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1586 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
1587
1588 angle = (b_A_size + b_B_size) / b_dist
1589
1590 delta_A = np.linalg.norm(np.array(s_A) - np.array(s_Ap))
1591
1592 if do_print:
1593 print("b_A_size =", b_A_size)
1594 print("b_B_size =", b_B_size)
1595 print("b_dist =", b_dist)
1596 print("angle =", angle)
1597
1598 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 ;) ).
1605 plt.figure()
1606 for order in range(1, 6):
1607 print("--------------------------------")
1608 print(f"Running FMM order = {order}")
1609 print("--------------------------------")
1610
1611 # set seed
1612 rng = np.random.default_rng(111)
1613
1614 N = 50000
1615
1616 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
1617 x_i_all = []
1618 for i in range(N):
1619 x_i_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1620
1621 # same for x_j
1622 x_j_all = []
1623 for i in range(N):
1624 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1625
1626 box_scale_fact_all = np.linspace(0, 0.1, N).tolist()
1627
1628 # same for box_1_center
1629 s_A_all = []
1630 s_Ap_all = []
1631 for p, box_scale_fact in zip(x_i_all, box_scale_fact_all):
1632 s_A_all.append(
1633 (
1634 p[0] + box_scale_fact * rng.uniform(-1, 1),
1635 p[1] + box_scale_fact * rng.uniform(-1, 1),
1636 p[2] + box_scale_fact * rng.uniform(-1, 1),
1637 )
1638 )
1639 s_Ap_all.append(
1640 (
1641 p[0] + box_scale_fact * rng.uniform(-1, 1),
1642 p[1] + box_scale_fact * rng.uniform(-1, 1),
1643 p[2] + box_scale_fact * rng.uniform(-1, 1),
1644 )
1645 )
1646
1647 # same for box_2_center
1648 s_B_all = []
1649 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1650 s_B_all.append(
1651 (
1652 p[0] + box_scale_fact * rng.uniform(-1, 1),
1653 p[1] + box_scale_fact * rng.uniform(-1, 1),
1654 p[2] + box_scale_fact * rng.uniform(-1, 1),
1655 )
1656 )
1657
1658 angles = []
1659 delta_A_all = []
1660 rel_errors = []
1661
1662 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):
1663 force_i, force_i_offset, angle, delta_A = test_grav_moment_offset(
1664 x_i, x_j, s_A, s_Ap, s_B, m_j, order, do_print=False
1665 )
1666
1667 abs_error = np.linalg.norm(force_i_offset - force_i)
1668 rel_error = abs_error / np.linalg.norm(force_i)
1669
1670 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
1671 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1672 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
1673 angle = (b_A_size + b_B_size) / b_dist
1674
1675 if angle > 5.0 or angle < 1e-4:
1676 continue
1677
1678 angles.append(angle)
1679 delta_A_all.append(delta_A)
1680 rel_errors.append(rel_error)
1681
1682 print(f"Computed for {len(angles)} cases")
1683
1684 plt.scatter(angles, rel_errors, s=1, label=f"FMM order = {order}")
1685
1686
1687 plt.xlabel("Angle")
1688 plt.ylabel("$|f_{\\rm fmm} - f_{\\rm fmm, offset}| / |f_{\\rm fmm}|$ (Relative error)")
1689 plt.xscale("log")
1690 plt.yscale("log")
1691 plt.title("Grav moment translation error")
1692 plt.legend(loc="lower right")
1693 plt.grid()
1694 plt.show()
--------------------------------
Running FMM order = 1
--------------------------------
Computed for 49965 cases
--------------------------------
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)
)
1760 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1761 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.519391031050157, v_001=0.870224438261286, 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.818204247354046, v_0001=3.478595136878889, 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.815100928798074, 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)
1771 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
1772 print("dM_k =", dM_k)
1773 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec((0, 0, 0))
1774 print("a_k =", a_k)
1775 force_i_fmm_sA_null = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
1776 Gconst = 1 # let's just set the grav constant to 1
1777 force_i_fmm_sA_null = -Gconst * np.array(force_i_fmm_sA_null)
1778 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.025579931908720904, v_02=0, v_11=1.012081300493465, v_12=0, v_22=0.9979099870995988),
t3=SymTensor3d_3(v_000=5.789190446027138, v_001=0.4640460581772785, v_002=0, v_011=-2.9347687627868018, v_012=0, v_022=-2.8544216832403366, v_111=-0.35343824590536155, v_112=0, v_122=-0.11060781227191739, v_222=0),
t4=SymTensor3d_4(v_0000=-17.09412455532619, 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.815100928798074, 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
1783 def analytic_force_i(x_i, x_j, Gconst):
1784 force_i_direct = (x_j[0] - x_i[0], x_j[1] - x_i[1], x_j[2] - x_i[2])
1785 force_i_direct /= np.linalg.norm(force_i_direct) ** 3
1786 force_i_direct *= m_i
1787 return force_i_direct
1788
1789
1790 force_i_exact = analytic_force_i(x_i, x_j, Gconst)
1791 print("force_i_exact =", force_i_exact)
force_i_exact = [-1. 0. 0.]
1794 abs_error = np.linalg.norm(force_i - force_i_exact)
1795 rel_error = abs_error / np.linalg.norm(force_i)
1796
1797
1798 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1799 b_dist = np.linalg.norm(np.array(x_i) - np.array(s_B))
1800 angle = (b_B_size) / b_dist
1801
1802 print("abs_error =", abs_error)
1803 print("rel_error =", rel_error)
1804 print("angle =", angle)
1805
1806 assert rel_error < 1e-2
abs_error = 0.0014495314137263013
rel_error = 0.001447602143490707
angle = 0.23249527748763862
Let’s code MM in a box
1812 def run_mm(x_i, x_j, s_B, m_j, order, do_print):
1813 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1814 r = (s_B[0] - x_i[0], s_B[1] - x_i[1], s_B[2] - x_i[2])
1815
1816 if do_print:
1817 print("x_i =", x_i)
1818 print("x_j =", x_j)
1819 print("s_B =", s_B)
1820 print("b_j =", b_j)
1821 print("r =", r)
1822
1823 # compute the tensor product of the displacment
1824 if order == 1:
1825 Q_n_B = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_j)
1826 elif order == 2:
1827 Q_n_B = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_j)
1828 elif order == 3:
1829 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
1830 elif order == 4:
1831 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
1832 elif order == 5:
1833 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
1834 else:
1835 raise ValueError("Invalid order")
1836
1837 # multiply by mass to get the mass moment
1838 Q_n_B *= m_j
1839
1840 if do_print:
1841 print("Q_n_B =", Q_n_B)
1842
1843 # green function gradients
1844 if order == 1:
1845 D_n = shamrock.phys.green_func_grav_cartesian_1_1(r)
1846 elif order == 2:
1847 D_n = shamrock.phys.green_func_grav_cartesian_1_2(r)
1848 elif order == 3:
1849 D_n = shamrock.phys.green_func_grav_cartesian_1_3(r)
1850 elif order == 4:
1851 D_n = shamrock.phys.green_func_grav_cartesian_1_4(r)
1852 elif order == 5:
1853 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1854 else:
1855 raise ValueError("Invalid order")
1856
1857 if do_print:
1858 print("D_n =", D_n)
1859
1860 if order == 1:
1861 result = shamrock.phys.contract_grav_moment_to_force_1(Q_n_B, D_n)
1862 elif order == 2:
1863 result = shamrock.phys.contract_grav_moment_to_force_2(Q_n_B, D_n)
1864 elif order == 3:
1865 result = shamrock.phys.contract_grav_moment_to_force_3(Q_n_B, D_n)
1866 elif order == 4:
1867 result = shamrock.phys.contract_grav_moment_to_force_4(Q_n_B, D_n)
1868 elif order == 5:
1869 result = shamrock.phys.contract_grav_moment_to_force_5(Q_n_B, D_n)
1870 else:
1871 raise ValueError("Invalid order")
1872
1873 Gconst = 1 # let's just set the grav constant to 1
1874 force_i = -Gconst * np.array(result)
1875 if do_print:
1876 print("force_i =", force_i)
1877
1878 force_i_exact = analytic_force_i(x_i, x_j, Gconst)
1879 if do_print:
1880 print("force_i_exact =", force_i_exact)
1881
1882 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1883 b_dist = np.linalg.norm(np.array(x_i) - np.array(s_B))
1884
1885 angle = (b_B_size) / b_dist
1886
1887 if do_print:
1888 print("b_A_size =", b_A_size)
1889 print("b_B_size =", b_B_size)
1890 print("b_dist =", b_dist)
1891 print("angle =", angle)
1892
1893 return force_i, force_i_exact, angle
1898 plt.figure()
1899 for order in range(1, 6):
1900 print("--------------------------------")
1901 print(f"Running MM order = {order}")
1902 print("--------------------------------")
1903
1904 # set seed
1905 rng = np.random.default_rng(111)
1906
1907 N = 50000
1908
1909 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
1910 x_i_all = []
1911 for i in range(N):
1912 x_i_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1913
1914 # same for x_j
1915 x_j_all = []
1916 for i in range(N):
1917 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1918
1919 box_scale_fact_all = np.linspace(0, 0.1, N).tolist()
1920
1921 # same for box_2_center
1922 s_B_all = []
1923 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1924 s_B_all.append(
1925 (
1926 p[0] + box_scale_fact * rng.uniform(-1, 1),
1927 p[1] + box_scale_fact * rng.uniform(-1, 1),
1928 p[2] + box_scale_fact * rng.uniform(-1, 1),
1929 )
1930 )
1931
1932 angles = []
1933 rel_errors = []
1934
1935 for x_i, x_j, s_B in zip(x_i_all, x_j_all, s_B_all):
1936 force_i, force_i_exact, angle = run_mm(x_i, x_j, s_B, m_j, order, do_print=False)
1937
1938 abs_error = np.linalg.norm(force_i - force_i_exact)
1939 rel_error = abs_error / np.linalg.norm(force_i_exact)
1940
1941 if angle > 5.0 or angle < 1e-4:
1942 continue
1943
1944 angles.append(angle)
1945 rel_errors.append(rel_error)
1946
1947 print(f"Computed for {len(angles)} cases")
1948
1949 plt.scatter(angles, rel_errors, s=1, label=f"MM order = {order}")
1950
1951
1952 def plot_powerlaw(order, center_y):
1953 X = [1e-3, 1e-2 / 3, 1e-1]
1954 Y = [center_y * (x) ** order for x in X]
1955 plt.plot(X, Y, linestyle="dashed", color="black")
1956 bbox = dict(boxstyle="round", fc="blanchedalmond", ec="orange", alpha=0.9)
1957 plt.text(X[1], Y[1], f"$\\propto x^{order}$", fontsize=9, bbox=bbox)
1958
1959
1960 plot_powerlaw(1, 1)
1961 plot_powerlaw(2, 1)
1962 plot_powerlaw(3, 1)
1963 plot_powerlaw(4, 1)
1964 plot_powerlaw(5, 1)
1965
1966 plt.xlabel("Angle")
1967 plt.ylabel("Relative Error")
1968 plt.xscale("log")
1969 plt.yscale("log")
1970 plt.title("MM precision")
1971 plt.legend(loc="lower right")
1972 plt.grid()
1973 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: (1 minutes 2.443 seconds)
Estimated memory usage: 202 MB