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
Use shamrock documentation style for matplotlib
19 shamrock.matplotlib.set_shamrock_mpl_style()
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
37 def draw_aabb(ax, aabb, color, alpha):
38 """
39 Draw a 3D AABB in matplotlib
40
41 Parameters
42 ----------
43 ax : matplotlib.Axes3D
44 The axis to draw the AABB on
45 aabb : shamrock.math.AABB_f64_3
46 The AABB to draw
47 color : str
48 The color of the AABB
49 alpha : float
50 The transparency of the AABB
51 """
52 xmin, ymin, zmin = aabb.lower
53 xmax, ymax, zmax = aabb.upper
54
55 points = [
56 aabb.lower,
57 (aabb.lower[0], aabb.lower[1], aabb.upper[2]),
58 (aabb.lower[0], aabb.upper[1], aabb.lower[2]),
59 (aabb.lower[0], aabb.upper[1], aabb.upper[2]),
60 (aabb.upper[0], aabb.lower[1], aabb.lower[2]),
61 (aabb.upper[0], aabb.lower[1], aabb.upper[2]),
62 (aabb.upper[0], aabb.upper[1], aabb.lower[2]),
63 aabb.upper,
64 ]
65
66 faces = [
67 [points[0], points[1], points[3], points[2]],
68 [points[4], points[5], points[7], points[6]],
69 [points[0], points[1], points[5], points[4]],
70 [points[2], points[3], points[7], points[6]],
71 [points[0], points[2], points[6], points[4]],
72 [points[1], points[3], points[7], points[5]],
73 ]
74
75 edges = [
76 [points[0], points[1]],
77 [points[0], points[2]],
78 [points[0], points[4]],
79 [points[1], points[3]],
80 [points[1], points[5]],
81 [points[2], points[3]],
82 [points[2], points[6]],
83 [points[3], points[7]],
84 [points[4], points[5]],
85 [points[4], points[6]],
86 [points[5], points[7]],
87 [points[6], points[7]],
88 ]
89
90 collection = Poly3DCollection(faces, alpha=alpha, color=color)
91 ax.add_collection3d(collection)
92
93 edge_collection = Line3DCollection(edges, color="k", alpha=alpha)
94 ax.add_collection3d(edge_collection)
95
96
97 def draw_arrow(ax, p1, p2, color, label, arr_scale=0.1):
98 length = np.linalg.norm(np.array(p2) - np.array(p1))
99 arrow_length_ratio = arr_scale / length
100 ax.quiver(
101 p1[0],
102 p1[1],
103 p1[2],
104 p2[0] - p1[0],
105 p2[1] - p1[1],
106 p2[2] - p1[2],
107 color=color,
108 label=label,
109 arrow_length_ratio=arrow_length_ratio,
110 )
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):
209 def plot_mass_moment_case(s_B, box_B_size, x_j):
210 box_B = shamrock.math.AABB_f64_3(
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 s_B[0] + box_B_size / 2,
218 s_B[1] + box_B_size / 2,
219 s_B[2] + box_B_size / 2,
220 ),
221 )
222
223 fig = plt.figure()
224 ax = fig.add_subplot(111, projection="3d")
225 ax.minorticks_off()
226
227 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
228
229 ax.scatter(s_B[0], s_B[1], s_B[2], color="black", label="s_B")
230
231 ax.scatter(x_j[0], x_j[1], x_j[2], color="red", label="$x_j$")
232
233 draw_aabb(ax, box_B, "blue", 0.2)
234
235 center_view = (0.0, 0.0, 0.0)
236 view_size = 2.0
237 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
238 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
239 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
240 ax.set_xlabel("X")
241 ax.set_ylabel("Y")
242 ax.set_zlabel("Z")
243 return ax
Let’s start with the following
254 s_B = (0, 0, 0)
255 box_B_size = 1
256
257 x_j = (0.2, 0.2, 0.2)
258 m_j = 1
259
260 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
261
262 ax = plot_mass_moment_case(s_B, box_B_size, x_j)
263 plt.title("Mass moment illustration")
264 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
265 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
294 s_B = (0, 0, 0)
295 box_B_size = 1
296
297 x_j = (0.5, 0.0, 0.0)
298 m_j = 1
299
300 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
301
302 ax = plot_mass_moment_case(s_B, box_B_size, x_j)
303 plt.title("Mass moment illustration")
304 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
305 plt.show()
306
307 Q_n_B = shamrock.math.SymTensorCollection_f64_0_5.from_vec(b_j)
308 Q_n_B *= m_j
309
310 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):
325 def plot_grav_moment_case(s_A, box_A_size, s_B, box_B_size, x_j):
326 box_A = shamrock.math.AABB_f64_3(
327 (
328 s_A[0] - box_A_size / 2,
329 s_A[1] - box_A_size / 2,
330 s_A[2] - box_A_size / 2,
331 ),
332 (
333 s_A[0] + box_A_size / 2,
334 s_A[1] + box_A_size / 2,
335 s_A[2] + box_A_size / 2,
336 ),
337 )
338
339 box_B = shamrock.math.AABB_f64_3(
340 (
341 s_B[0] - box_B_size / 2,
342 s_B[1] - box_B_size / 2,
343 s_B[2] - box_B_size / 2,
344 ),
345 (
346 s_B[0] + box_B_size / 2,
347 s_B[1] + box_B_size / 2,
348 s_B[2] + box_B_size / 2,
349 ),
350 )
351
352 fig = plt.figure()
353 ax = fig.add_subplot(111, projection="3d")
354 ax.minorticks_off()
355
356 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
357
358 draw_arrow(ax, s_A, s_B, "purple", "$r = s_B - s_A$")
359
360 ax.scatter(s_A[0], s_A[1], s_A[2], color="black", label="s_A")
361
362 ax.scatter(s_B[0], s_B[1], s_B[2], color="green", label="s_B")
363
364 ax.scatter(x_j[0], x_j[1], x_j[2], color="red", label="$x_j$")
365
366 draw_aabb(ax, box_A, "blue", 0.1)
367 draw_aabb(ax, box_B, "red", 0.1)
368
369 center_view = (0.5, 0.0, 0.0)
370 view_size = 2.0
371 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
372 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
373 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
374 ax.set_xlabel("X")
375 ax.set_ylabel("Y")
376 ax.set_zlabel("Z")
377 return ax
Let’s now show the example of a gravitational moment, for the following case
387 s_B = (0, 0, 0)
388 s_A = (1, 0, 0)
389
390 box_B_size = 0.5
391 box_A_size = 0.5
392
393 x_j = (0.2, 0.2, 0.0)
394 m_j = 1
395
396 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
397 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
398
399 ax = plot_grav_moment_case(s_A, box_A_size, s_B, box_B_size, x_j)
400 plt.title("Grav moment illustration")
401 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
402 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
412 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
413 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
417 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
418 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):
432 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):
433 box_A = shamrock.math.AABB_f64_3(
434 (
435 s_A[0] - box_A_size / 2,
436 s_A[1] - box_A_size / 2,
437 s_A[2] - box_A_size / 2,
438 ),
439 (
440 s_A[0] + box_A_size / 2,
441 s_A[1] + box_A_size / 2,
442 s_A[2] + box_A_size / 2,
443 ),
444 )
445
446 box_B = shamrock.math.AABB_f64_3(
447 (
448 s_B[0] - box_B_size / 2,
449 s_B[1] - box_B_size / 2,
450 s_B[2] - box_B_size / 2,
451 ),
452 (
453 s_B[0] + box_B_size / 2,
454 s_B[1] + box_B_size / 2,
455 s_B[2] + box_B_size / 2,
456 ),
457 )
458
459 fig = plt.figure()
460 ax = fig.add_subplot(111, projection="3d")
461 ax.minorticks_off()
462
463 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
464 draw_arrow(ax, s_A, x_i, "blue", "$a_i = x_i - s_A$")
465
466 draw_arrow(ax, s_A, s_B, "purple", "$r = s_B - s_A$")
467
468 ax.scatter(s_A[0], s_A[1], s_A[2], color="black", label="s_A")
469
470 ax.scatter(s_B[0], s_B[1], s_B[2], color="green", label="s_B")
471
472 ax.scatter(x_i[0], x_i[1], x_i[2], color="orange", label="$x_i$")
473
474 ax.scatter(x_j[0], x_j[1], x_j[2], color="red", label="$x_j$")
475
476 draw_arrow(ax, x_i, x_i + force_i * fscale_fact, "green", "$f_i$")
477 draw_arrow(ax, x_i, x_i + force_i_exact * fscale_fact, "red", "$f_i$ (exact)")
478
479 abs_error = np.linalg.norm(force_i - force_i_exact)
480 rel_error = abs_error / np.linalg.norm(force_i_exact)
481
482 draw_aabb(ax, box_A, "blue", 0.1)
483 draw_aabb(ax, box_B, "red", 0.1)
484
485 center_view = (0.5, 0.0, 0.0)
486 view_size = 2.0
487 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
488 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
489 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
490 ax.set_xlabel("X")
491 ax.set_ylabel("Y")
492 ax.set_zlabel("Z")
493
494 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)
507 s_B = (0, 0, 0)
508 s_A = (1, 0, 0)
509
510 box_B_size = 0.5
511 box_A_size = 0.5
512
513 x_j = (0.2, 0.2, 0.0)
514 x_i = (1.2, 0.2, 0.0)
515 m_j = 1
516 m_i = 1
517
518 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
519 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
520 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)
)
528 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
529 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)
)
532 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
533 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)
)
536 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_i)
537 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
548 def analytic_force_i(x_i, x_j, Gconst):
549 force_i_direct = (x_j[0] - x_i[0], x_j[1] - x_i[1], x_j[2] - x_i[2])
550 force_i_direct /= np.linalg.norm(force_i_direct) ** 3
551 force_i_direct *= m_i
552 return force_i_direct
553
554
555 force_i_exact = analytic_force_i(x_i, x_j, Gconst)
556 print("force_i_exact =", force_i_exact)
force_i_exact = [-1. 0. 0.]
This yields the following case
560 ax, rel_error, abs_error = plot_fmm_case(
561 s_A, box_A_size, x_i, s_B, box_B_size, x_j, force_i, force_i_exact, 0.5
562 )
563
564 plt.title(f"FMM, rel error={rel_error}")
565 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
566 plt.show()
567
568 print("force_i =", force_i)
569 print("force_i_exact =", force_i_exact)
570 print("abs error =", abs_error)
571 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
584 def run_fmm(x_i, x_j, s_A, s_B, m_j, order, do_print):
585 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
586 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
587 a_i = (x_i[0] - s_A[0], x_i[1] - s_A[1], x_i[2] - s_A[2])
588
589 if do_print:
590 print("x_i =", x_i)
591 print("x_j =", x_j)
592 print("s_A =", s_A)
593 print("s_B =", s_B)
594 print("b_j =", b_j)
595 print("r =", r)
596 print("a_i =", a_i)
597
598 # compute the tensor product of the displacment
599 if order == 1:
600 Q_n_B = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_j)
601 elif order == 2:
602 Q_n_B = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_j)
603 elif order == 3:
604 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
605 elif order == 4:
606 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
607 elif order == 5:
608 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
609 else:
610 raise ValueError("Invalid order")
611
612 # multiply by mass to get the mass moment
613 Q_n_B *= m_j
614
615 if do_print:
616 print("Q_n_B =", Q_n_B)
617
618 # green function gradients
619 if order == 1:
620 D_n = shamrock.phys.green_func_grav_cartesian_1_1(r)
621 elif order == 2:
622 D_n = shamrock.phys.green_func_grav_cartesian_1_2(r)
623 elif order == 3:
624 D_n = shamrock.phys.green_func_grav_cartesian_1_3(r)
625 elif order == 4:
626 D_n = shamrock.phys.green_func_grav_cartesian_1_4(r)
627 elif order == 5:
628 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
629 else:
630 raise ValueError("Invalid order")
631
632 if do_print:
633 print("D_n =", D_n)
634
635 if order == 1:
636 dM_k = shamrock.phys.get_dM_mat_1(D_n, Q_n_B)
637 elif order == 2:
638 dM_k = shamrock.phys.get_dM_mat_2(D_n, Q_n_B)
639 elif order == 3:
640 dM_k = shamrock.phys.get_dM_mat_3(D_n, Q_n_B)
641 elif order == 4:
642 dM_k = shamrock.phys.get_dM_mat_4(D_n, Q_n_B)
643 elif order == 5:
644 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
645 else:
646 raise ValueError("Invalid order")
647
648 if do_print:
649 print("dM_k =", dM_k)
650
651 if order == 1:
652 a_k = shamrock.math.SymTensorCollection_f64_0_0.from_vec(a_i)
653 elif order == 2:
654 a_k = shamrock.math.SymTensorCollection_f64_0_1.from_vec(a_i)
655 elif order == 3:
656 a_k = shamrock.math.SymTensorCollection_f64_0_2.from_vec(a_i)
657 elif order == 4:
658 a_k = shamrock.math.SymTensorCollection_f64_0_3.from_vec(a_i)
659 elif order == 5:
660 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_i)
661 else:
662 raise ValueError("Invalid order")
663
664 if do_print:
665 print("a_k =", a_k)
666
667 if order == 1:
668 result = shamrock.phys.contract_grav_moment_to_force_1(a_k, dM_k)
669 elif order == 2:
670 result = shamrock.phys.contract_grav_moment_to_force_2(a_k, dM_k)
671 elif order == 3:
672 result = shamrock.phys.contract_grav_moment_to_force_3(a_k, dM_k)
673 elif order == 4:
674 result = shamrock.phys.contract_grav_moment_to_force_4(a_k, dM_k)
675 elif order == 5:
676 result = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
677 else:
678 raise ValueError("Invalid order")
679
680 Gconst = 1 # let's just set the grav constant to 1
681 force_i = -Gconst * np.array(result)
682 if do_print:
683 print("force_i =", force_i)
684
685 force_i_exact = analytic_force_i(x_i, x_j, Gconst)
686 if do_print:
687 print("force_i_exact =", force_i_exact)
688
689 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
690 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
691 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
692
693 angle = (b_A_size + b_B_size) / b_dist
694
695 if do_print:
696 print("b_A_size =", b_A_size)
697 print("b_B_size =", b_B_size)
698 print("b_dist =", b_dist)
699 print("angle =", angle)
700
701 return force_i, force_i_exact, angle
Let’s try with some new parameters
706 s_B = (0, 0, 0)
707 s_A = (1, 0, 0)
708
709 box_B_size = 0.5
710 box_A_size = 0.5
711
712 x_j = (0.2, 0.2, 0.0)
713 x_i = (1.2, 0.2, 0.2)
714 m_j = 1
715 m_i = 1
716
717 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order=5, do_print=True)
718 ax, rel_error, abs_error = plot_fmm_case(
719 s_A, box_A_size, x_i, s_B, box_B_size, x_j, force_i, force_i_exact, 0.5
720 )
721
722 plt.title(f"FMM angle={angle:.5f} rel error={rel_error:.2e}")
723 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
724 plt.show()
725
726 print("force_i =", force_i)
727 print("force_i_exact =", force_i_exact)
728 print("abs error =", abs_error)
729 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#
736 s_B = (0, 0, 0)
737 s_A = (1, 0, 0)
738
739 box_B_size = 0.5
740 box_A_size = 0.5
741
742 x_j = (0.2, 0.2, 0.0)
743 x_i = (0.8, 0.2, 0.2)
744 m_j = 1
745 m_i = 1
746
747
748 for order in range(1, 6):
749 print("--------------------------------")
750 print(f"Running FMM order = {order}")
751 print("--------------------------------")
752
753 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order, do_print=True)
754 ax, rel_error, abs_error = plot_fmm_case(
755 s_A, box_A_size, x_i, s_B, box_B_size, x_j, force_i, force_i_exact, 0.2
756 )
757
758 plt.title(f"FMM order={order} angle={angle:.5f} rel error={rel_error:.2e}")
759 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
760 plt.show()
761
762 print("force_i =", force_i)
763 print("force_i_exact =", force_i_exact)
764 print("abs error =", abs_error)
765 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.5832193498525184
rel error = 0.6332877399410074
--------------------------------
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.8219626217344302
rel error = 0.3287850486937721
--------------------------------
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.38198731183411494
rel error = 0.15279492473364598
--------------------------------
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.17077083634710036
rel error = 0.06830833453884014
--------------------------------
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.07940820523782473
rel error = 0.031763282095129894
Sweeping through angles#
773 s_B = (0, 0, 0)
774 s_A_all = [(0.8, 0, 0), (1, 0, 0), (1.2, 0, 0)]
775
776 box_B_size = 0.5
777 box_A_size = 0.5
778
779 x_j = (0.2, 0.2, 0.0)
780 x_i = (0.8, 0.2, 0.2)
781 m_j = 1
782 m_i = 1
783
784 order = 3
785
786 for s_A in s_A_all:
787 print("--------------------------------")
788 print(f"Running FMM s_a = {s_A}")
789 print("--------------------------------")
790
791 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order, do_print=True)
792 ax, rel_error, abs_error = plot_fmm_case(
793 s_A, box_A_size, x_i, s_B, box_B_size, x_j, force_i, force_i_exact, 0.2
794 )
795
796 plt.title(f"FMM order={order} angle={angle:.5f} rel error={rel_error:.2e}")
797 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
798 plt.show()
799
800 print("force_i =", force_i)
801 print("force_i_exact =", force_i_exact)
802 print("abs error =", abs_error)
803 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.15966288352037034
rel error = 0.06386515340814813
--------------------------------
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.38198731183411494
rel error = 0.15279492473364598
--------------------------------
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.7015858309588155
rel error = 0.2806343323835262
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.
817 plt.figure()
818 for order in range(1, 6):
819 print("--------------------------------")
820 print(f"Running FMM order = {order}")
821 print("--------------------------------")
822
823 # set seed
824 rng = np.random.default_rng(111)
825
826 N = 50000
827
828 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
829 x_i_all = []
830 for i in range(N):
831 x_i_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
832
833 # same for x_j
834 x_j_all = []
835 for i in range(N):
836 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
837
838 box_scale_fact_all = np.linspace(0, 0.1, N).tolist()
839
840 # same for box_1_center
841 s_A_all = []
842 for p, box_scale_fact in zip(x_i_all, box_scale_fact_all):
843 s_A_all.append(
844 (
845 p[0] + box_scale_fact * rng.uniform(-1, 1),
846 p[1] + box_scale_fact * rng.uniform(-1, 1),
847 p[2] + box_scale_fact * rng.uniform(-1, 1),
848 )
849 )
850
851 # same for box_2_center
852 s_B_all = []
853 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
854 s_B_all.append(
855 (
856 p[0] + box_scale_fact * rng.uniform(-1, 1),
857 p[1] + box_scale_fact * rng.uniform(-1, 1),
858 p[2] + box_scale_fact * rng.uniform(-1, 1),
859 )
860 )
861
862 angles = []
863 rel_errors = []
864
865 for x_i, x_j, s_A, s_B in zip(x_i_all, x_j_all, s_A_all, s_B_all):
866 force_i, force_i_exact, angle = run_fmm(x_i, x_j, s_A, s_B, m_j, order, do_print=False)
867
868 abs_error = np.linalg.norm(force_i - force_i_exact)
869 rel_error = abs_error / np.linalg.norm(force_i_exact)
870
871 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
872 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
873 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
874 angle = (b_A_size + b_B_size) / b_dist
875
876 if angle > 5.0 or angle < 1e-4:
877 continue
878
879 angles.append(angle)
880 rel_errors.append(rel_error)
881
882 print(f"Computed for {len(angles)} cases")
883
884 plt.scatter(angles, rel_errors, s=1, label=f"FMM order = {order}")
885
886
887 def plot_powerlaw(order, center_y):
888 X = [1e-3, 1e-2 / 3, 1e-1]
889 Y = [center_y * (x) ** order for x in X]
890 plt.plot(X, Y, linestyle="dashed", color="black")
891 bbox = dict(boxstyle="round", fc="blanchedalmond", ec="orange", alpha=0.9)
892 plt.text(X[1], Y[1], f"$\\propto x^{order}$", fontsize=9, bbox=bbox)
893
894
895 plot_powerlaw(1, 1)
896 plot_powerlaw(2, 1)
897 plot_powerlaw(3, 1)
898 plot_powerlaw(4, 1)
899 plot_powerlaw(5, 1)
900
901 plt.xlabel("Angle")
902 plt.ylabel("Relative Error")
903 plt.xscale("log")
904 plt.yscale("log")
905 plt.title("FMM precision")
906 plt.legend(loc="lower right")
907 plt.grid()
908 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. Therefore 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.
938 s_B = (0, 0, 0)
939 box_B_size = 1.0
940 x_j = (0.2, 0.2, 0.0)
941 m_j = 1
942
943 s_B_new = (0.3, 0.3, 0.3)
def plot_mass_moment_offset(s_B, s_B_new, box_B_size):
952 def plot_mass_moment_offset(s_B, s_B_new, box_B_size):
953 box_B = shamrock.math.AABB_f64_3(
954 (
955 s_B[0] - box_B_size / 2,
956 s_B[1] - box_B_size / 2,
957 s_B[2] - box_B_size / 2,
958 ),
959 (
960 s_B[0] + box_B_size / 2,
961 s_B[1] + box_B_size / 2,
962 s_B[2] + box_B_size / 2,
963 ),
964 )
965
966 box_B_new = shamrock.math.AABB_f64_3(
967 (
968 s_B_new[0] - box_B_size / 2,
969 s_B_new[1] - box_B_size / 2,
970 s_B_new[2] - box_B_size / 2,
971 ),
972 (
973 s_B_new[0] + box_B_size / 2,
974 s_B_new[1] + box_B_size / 2,
975 s_B_new[2] + box_B_size / 2,
976 ),
977 )
978
979 fig = plt.figure()
980 ax = fig.add_subplot(111, projection="3d")
981 ax.minorticks_off()
982
983 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
984 draw_arrow(ax, s_B_new, x_j, "red", "$b_j' = x_j - s_B'$")
985
986 ax.scatter(s_B[0], s_B[1], s_B[2], color="black", label="s_B")
987 ax.scatter(s_B_new[0], s_B_new[1], s_B_new[2], color="red", label="s_B'")
988 ax.scatter(x_j[0], x_j[1], x_j[2], color="blue", label="$x_j$")
989
990 draw_aabb(ax, box_B, "blue", 0.2)
991 draw_aabb(ax, box_B_new, "red", 0.2)
992
993 center_view = (0.0, 0.0, 0.0)
994 view_size = 2.0
995 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
996 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
997 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
998 ax.set_xlabel("X")
999 ax.set_ylabel("Y")
1000 ax.set_zlabel("Z")
1001
1002 return ax
1012 plot_mass_moment_offset(s_B, s_B_new, box_B_size)
1013
1014 plt.title("Mass moment offset illustration")
1015 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
1016 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’
1036 Q_n_B_offset = shamrock.phys.offset_multipole_5(Q_n_B, s_B, s_B_new)
1037 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’
1043 def tensor_collect_norm(d):
1044 # detect the type of the tensor collection
1045 if isinstance(d, shamrock.math.SymTensorCollection_f64_0_5):
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 + np.sqrt(d.t5.inner(d.t5)) / 120
1053 )
1054 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_4):
1055 return (
1056 np.sqrt(d.t0 * d.t0)
1057 + np.sqrt(d.t1.inner(d.t1))
1058 + np.sqrt(d.t2.inner(d.t2)) / 2
1059 + np.sqrt(d.t3.inner(d.t3)) / 6
1060 + np.sqrt(d.t4.inner(d.t4)) / 24
1061 )
1062 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_3):
1063 return (
1064 np.sqrt(d.t0 * d.t0)
1065 + np.sqrt(d.t1.inner(d.t1))
1066 + np.sqrt(d.t2.inner(d.t2)) / 2
1067 + np.sqrt(d.t3.inner(d.t3)) / 6
1068 )
1069 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_2):
1070 return np.sqrt(d.t0 * d.t0) + np.sqrt(d.t1.inner(d.t1)) + np.sqrt(d.t2.inner(d.t2)) / 2
1071 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_1):
1072 return np.sqrt(d.t0 * d.t0) + np.sqrt(d.t1.inner(d.t1))
1073 elif isinstance(d, shamrock.math.SymTensorCollection_f64_0_0):
1074 return np.sqrt(d.t0 * d.t0)
1075 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_5):
1076 return (
1077 np.sqrt(d.t1.inner(d.t1))
1078 + np.sqrt(d.t2.inner(d.t2)) / 2
1079 + np.sqrt(d.t3.inner(d.t3)) / 6
1080 + np.sqrt(d.t4.inner(d.t4)) / 24
1081 + np.sqrt(d.t5.inner(d.t5)) / 120
1082 )
1083 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_4):
1084 return (
1085 np.sqrt(d.t1.inner(d.t1))
1086 + np.sqrt(d.t2.inner(d.t2)) / 2
1087 + np.sqrt(d.t3.inner(d.t3)) / 6
1088 + np.sqrt(d.t4.inner(d.t4)) / 24
1089 )
1090 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_3):
1091 return (
1092 np.sqrt(d.t1.inner(d.t1))
1093 + np.sqrt(d.t2.inner(d.t2)) / 2
1094 + np.sqrt(d.t3.inner(d.t3)) / 6
1095 )
1096 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_2):
1097 return np.sqrt(d.t1.inner(d.t1)) + np.sqrt(d.t2.inner(d.t2)) / 2
1098 elif isinstance(d, shamrock.math.SymTensorCollection_f64_1_1):
1099 return np.sqrt(d.t1.inner(d.t1))
1100 else:
1101 raise ValueError(f"Unsupported tensor collection type: {type(d)}")
1102
1103
1104 print("Q_n_B norm =", tensor_collect_norm(Q_n_B))
1105 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
1109 delta = Q_n_B_offset - Q_n_Bp
1110 print("delta =", delta)
1111
1112
1113 sqdist_t0 = delta.t0 * delta.t0
1114 sqdist_t1 = delta.t1.inner(delta.t1)
1115 sqdist_t2 = delta.t2.inner(delta.t2)
1116 sqdist_t3 = delta.t3.inner(delta.t3)
1117 sqdist_t4 = delta.t4.inner(delta.t4)
1118 sqdist_t5 = delta.t5.inner(delta.t5)
1119 print("sqdist_t0 =", sqdist_t0)
1120 print("sqdist_t1 =", sqdist_t1)
1121 print("sqdist_t2 =", sqdist_t2)
1122 print("sqdist_t3 =", sqdist_t3)
1123 print("sqdist_t4 =", sqdist_t4)
1124 print("sqdist_t5 =", sqdist_t5)
1125
1126 norm_delta = (
1127 np.sqrt(sqdist_t0)
1128 + np.sqrt(sqdist_t1)
1129 + np.sqrt(sqdist_t2) / 2
1130 + np.sqrt(sqdist_t3) / 6
1131 + np.sqrt(sqdist_t4) / 24
1132 + np.sqrt(sqdist_t5) / 120
1133 )
1134 print("norm_delta =", norm_delta)
1135
1136 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
1141 plt.figure()
1142
1143 for order in range(0, 6):
1144 # set seed
1145 rng = np.random.default_rng(111)
1146
1147 N = 50000
1148
1149 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
1150 x_j_all = []
1151 for i in range(N):
1152 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1153
1154 box_scale_fact_all = np.linspace(0, 1, N).tolist()
1155
1156 # same for box_1_center
1157 s_B_all = []
1158 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1159 s_B_all.append(
1160 (
1161 p[0] + box_scale_fact * rng.uniform(-1, 1),
1162 p[1] + box_scale_fact * rng.uniform(-1, 1),
1163 p[2] + box_scale_fact * rng.uniform(-1, 1),
1164 )
1165 )
1166
1167 # same for box_2_center
1168 s_Bp_all = []
1169 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1170 s_Bp_all.append(
1171 (
1172 p[0] + box_scale_fact * rng.uniform(-1, 1),
1173 p[1] + box_scale_fact * rng.uniform(-1, 1),
1174 p[2] + box_scale_fact * rng.uniform(-1, 1),
1175 )
1176 )
1177
1178 center_distances = []
1179 rel_errors = []
1180 for x_j, s_B, s_Bp in zip(x_j_all, s_B_all, s_Bp_all):
1181 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1182 b_jp = (x_j[0] - s_Bp[0], x_j[1] - s_Bp[1], x_j[2] - s_Bp[2])
1183
1184 if order == 5:
1185 Q_n_B = shamrock.math.SymTensorCollection_f64_0_5.from_vec(b_j)
1186 Q_n_B *= m_j
1187
1188 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_5.from_vec(b_jp)
1189 Q_n_Bp *= m_j
1190
1191 Q_n_B_offset = shamrock.phys.offset_multipole_5(Q_n_B, s_B, s_Bp)
1192 elif order == 4:
1193 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
1194 Q_n_B *= m_j
1195
1196 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_jp)
1197 Q_n_Bp *= m_j
1198
1199 Q_n_B_offset = shamrock.phys.offset_multipole_4(Q_n_B, s_B, s_Bp)
1200 elif order == 3:
1201 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
1202 Q_n_B *= m_j
1203
1204 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_jp)
1205 Q_n_Bp *= m_j
1206
1207 Q_n_B_offset = shamrock.phys.offset_multipole_3(Q_n_B, s_B, s_Bp)
1208 elif order == 2:
1209 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
1210 Q_n_B *= m_j
1211
1212 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_jp)
1213 Q_n_Bp *= m_j
1214
1215 Q_n_B_offset = shamrock.phys.offset_multipole_2(Q_n_B, s_B, s_Bp)
1216 elif order == 1:
1217 Q_n_B = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_j)
1218 Q_n_B *= m_j
1219
1220 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_jp)
1221 Q_n_Bp *= m_j
1222
1223 Q_n_B_offset = shamrock.phys.offset_multipole_1(Q_n_B, s_B, s_Bp)
1224 elif order == 0:
1225 Q_n_B = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_j)
1226 Q_n_B *= m_j
1227
1228 Q_n_Bp = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_jp)
1229 Q_n_Bp *= m_j
1230
1231 Q_n_B_offset = shamrock.phys.offset_multipole_0(Q_n_B, s_B, s_Bp)
1232 else:
1233 raise ValueError(f"Unsupported offset order: {order}")
1234
1235 delta = Q_n_B_offset - Q_n_Bp
1236
1237 rel_error = tensor_collect_norm(delta) / tensor_collect_norm(Q_n_B)
1238 rel_errors.append(rel_error)
1239
1240 center_distances.append(np.linalg.norm(np.array(s_B) - np.array(s_Bp)))
1241
1242 plt.scatter(center_distances, rel_errors, s=1, label=f"multipole order = {order}")
1243
1244 plt.xlabel("$\\vert \\vert s_B - s_B'\\vert \\vert$")
1245 plt.ylabel(
1246 "$\\vert \\vert Q_n(s_B) - Q_n(s_B') \\vert \\vert / \\vert \\vert Q_n(s_B) \\vert \\vert$ (relative error) "
1247 )
1248 plt.xscale("log")
1249 plt.yscale("log")
1250 plt.title("Mass moment offset precision")
1251 plt.legend(loc="lower right")
1252 plt.grid()
1253 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.
1282 def plot_grav_moment_offset(s_A, s_Ap, s_B, box_A_size, box_B_size, x_j):
1283 box_A = shamrock.math.AABB_f64_3(
1284 (
1285 s_A[0] - box_A_size / 2,
1286 s_A[1] - box_A_size / 2,
1287 s_A[2] - box_A_size / 2,
1288 ),
1289 (
1290 s_A[0] + box_A_size / 2,
1291 s_A[1] + box_A_size / 2,
1292 s_A[2] + box_A_size / 2,
1293 ),
1294 )
1295
1296 box_Ap = shamrock.math.AABB_f64_3(
1297 (
1298 s_Ap[0] - box_A_size / 2,
1299 s_Ap[1] - box_A_size / 2,
1300 s_Ap[2] - box_A_size / 2,
1301 ),
1302 (
1303 s_Ap[0] + box_A_size / 2,
1304 s_Ap[1] + box_A_size / 2,
1305 s_Ap[2] + box_A_size / 2,
1306 ),
1307 )
1308
1309 box_B = shamrock.math.AABB_f64_3(
1310 (
1311 s_B[0] - box_B_size / 2,
1312 s_B[1] - box_B_size / 2,
1313 s_B[2] - box_B_size / 2,
1314 ),
1315 (
1316 s_B[0] + box_B_size / 2,
1317 s_B[1] + box_B_size / 2,
1318 s_B[2] + box_B_size / 2,
1319 ),
1320 )
1321
1322 fig = plt.figure()
1323 ax = fig.add_subplot(111, projection="3d")
1324 ax.minorticks_off()
1325
1326 draw_arrow(ax, s_A, s_B, "purple", "$r = s_B - s_A$")
1327 draw_arrow(ax, s_Ap, s_B, "purple", "$r' = s_B - s_A'$")
1328
1329 ax.scatter(s_A[0], s_A[1], s_A[2], color="black", label="s_A")
1330 ax.scatter(s_Ap[0], s_Ap[1], s_Ap[2], color="black", label="s_Ap")
1331 ax.scatter(s_B[0], s_B[1], s_B[2], color="black", label="s_B")
1332
1333 draw_arrow(ax, s_B, x_j, "black", "$b_j = x_j - s_B$")
1334
1335 ax.scatter(x_j[0], x_j[1], x_j[2], color="red", label="$x_j$")
1336
1337 draw_aabb(ax, box_A, "blue", 0.1)
1338 draw_aabb(ax, box_Ap, "cyan", 0.1)
1339 draw_aabb(ax, box_B, "red", 0.1)
1340
1341 center_view = (0.5, 0.0, 0.0)
1342 view_size = 2.0
1343 ax.set_xlim(center_view[0] - view_size / 2, center_view[0] + view_size / 2)
1344 ax.set_ylim(center_view[1] - view_size / 2, center_view[1] + view_size / 2)
1345 ax.set_zlim(center_view[2] - view_size / 2, center_view[2] + view_size / 2)
1346 ax.set_xlabel("X")
1347 ax.set_ylabel("Y")
1348 ax.set_zlabel("Z")
1349
1350 return ax, rel_error, abs_error
1351
1352
1353 s_B = (0, -0.2, -0.2)
1354 s_A = (1, 0, 0)
1355 s_Ap = (1.1, 0.1, 0.0)
1356
1357 box_B_size = 0.5
1358 box_A_size = 0.5
1359
1360 x_j = (0.2, 0.0, -0.5)
1361 x_i = (1.2, 0.2, 0.0)
1362 m_j = 1
1363 m_i = 1
1364
1365 plot_grav_moment_offset(s_A, s_Ap, s_B, box_A_size, box_B_size, x_j)
1366
1367 plt.title("Mass moment offset illustration")
1368 plt.legend(loc="center left", bbox_to_anchor=(-0.3, 0.5))
1369 plt.show()
1370
1371 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1372 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
1373 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)
)
1382 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1383 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
1384 # print("D_n =",D_n)
1385 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)
)
1388 Dp_n = shamrock.phys.green_func_grav_cartesian_1_5(rp)
1389 dMp_k = shamrock.phys.get_dM_mat_5(Dp_n, Q_n_B)
1390 # print("Dp_n =",Dp_n)
1391 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
1395 dM_k_offset = shamrock.phys.offset_dM_mat_5(dM_k, s_A, s_Ap)
1396 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.
1403 delta = dM_k_offset - dMp_k
1404
1405 print("delta =", delta)
1406 print("sqdist_t1 =", np.sqrt(delta.t1.inner(delta.t1)))
1407 print("sqdist_t2 =", np.sqrt(delta.t2.inner(delta.t2)) / 2)
1408 print("sqdist_t3 =", np.sqrt(delta.t3.inner(delta.t3)) / 6)
1409 print("sqdist_t4 =", np.sqrt(delta.t4.inner(delta.t4)) / 24)
1410 print("sqdist_t5 =", np.sqrt(delta.t5.inner(delta.t5)) / 120)
1411 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
1415 a_i = (x_i[0] - s_A[0], x_i[1] - s_A[1], x_i[2] - s_A[2])
1416 a_ip = (x_i[0] - s_Ap[0], x_i[1] - s_Ap[1], x_i[2] - s_Ap[2])
1417
1418 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_i)
1419 a_kp = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_ip)
1420
1421 print("a_k =", a_k)
1422 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)
)
1425 result = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
1426 resultp = shamrock.phys.contract_grav_moment_to_force_5(a_kp, dMp_k)
1427 result_offset = shamrock.phys.contract_grav_moment_to_force_5(a_kp, dM_k_offset)
1428
1429 print("force_i =", -Gconst * np.array(result))
1430 print("force_ip =", -Gconst * np.array(resultp))
1431 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
1435 delta_f = np.linalg.norm(np.array(result_offset) - np.array(result))
1436 delta_f /= np.linalg.norm(np.array(result))
1437 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
1442 def test_grav_moment_offset(x_i, x_j, s_A, s_Ap, s_B, m_j, order, do_print):
1443 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1444 r = (s_B[0] - s_A[0], s_B[1] - s_A[1], s_B[2] - s_A[2])
1445 a_i = (x_i[0] - s_A[0], x_i[1] - s_A[1], x_i[2] - s_A[2])
1446 a_ip = (x_i[0] - s_Ap[0], x_i[1] - s_Ap[1], x_i[2] - s_Ap[2])
1447
1448 if do_print:
1449 print("x_i =", x_i)
1450 print("x_j =", x_j)
1451 print("s_A =", s_A)
1452 print("s_Ap =", s_Ap)
1453 print("s_B =", s_B)
1454 print("b_j =", b_j)
1455 print("r =", r)
1456 print("a_i =", a_i)
1457
1458 # compute the tensor product of the displacment
1459 if order == 1:
1460 Q_n_B = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_j)
1461 elif order == 2:
1462 Q_n_B = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_j)
1463 elif order == 3:
1464 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
1465 elif order == 4:
1466 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
1467 elif order == 5:
1468 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
1469 else:
1470 raise ValueError("Invalid order")
1471
1472 # multiply by mass to get the mass moment
1473 Q_n_B *= m_j
1474
1475 if do_print:
1476 print("Q_n_B =", Q_n_B)
1477
1478 # green function gradients
1479 if order == 1:
1480 D_n = shamrock.phys.green_func_grav_cartesian_1_1(r)
1481 elif order == 2:
1482 D_n = shamrock.phys.green_func_grav_cartesian_1_2(r)
1483 elif order == 3:
1484 D_n = shamrock.phys.green_func_grav_cartesian_1_3(r)
1485 elif order == 4:
1486 D_n = shamrock.phys.green_func_grav_cartesian_1_4(r)
1487 elif order == 5:
1488 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1489 else:
1490 raise ValueError("Invalid order")
1491
1492 if do_print:
1493 print("D_n =", D_n)
1494
1495 if order == 1:
1496 dM_k = shamrock.phys.get_dM_mat_1(D_n, Q_n_B)
1497 elif order == 2:
1498 dM_k = shamrock.phys.get_dM_mat_2(D_n, Q_n_B)
1499 elif order == 3:
1500 dM_k = shamrock.phys.get_dM_mat_3(D_n, Q_n_B)
1501 elif order == 4:
1502 dM_k = shamrock.phys.get_dM_mat_4(D_n, Q_n_B)
1503 elif order == 5:
1504 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
1505 else:
1506 raise ValueError("Invalid order")
1507
1508 if do_print:
1509 print("dM_k =", dM_k)
1510
1511 if order == 5:
1512 dM_k_offset = shamrock.phys.offset_dM_mat_5(dM_k, s_A, s_Ap)
1513 elif order == 4:
1514 dM_k_offset = shamrock.phys.offset_dM_mat_4(dM_k, s_A, s_Ap)
1515 elif order == 3:
1516 dM_k_offset = shamrock.phys.offset_dM_mat_3(dM_k, s_A, s_Ap)
1517 elif order == 2:
1518 dM_k_offset = shamrock.phys.offset_dM_mat_2(dM_k, s_A, s_Ap)
1519 elif order == 1:
1520 dM_k_offset = shamrock.phys.offset_dM_mat_1(dM_k, s_A, s_Ap)
1521 else:
1522 raise ValueError("Invalid order")
1523
1524 if do_print:
1525 print("dM_k_offset =", dM_k_offset)
1526
1527 if order == 1:
1528 a_k = shamrock.math.SymTensorCollection_f64_0_0.from_vec(a_i)
1529 elif order == 2:
1530 a_k = shamrock.math.SymTensorCollection_f64_0_1.from_vec(a_i)
1531 elif order == 3:
1532 a_k = shamrock.math.SymTensorCollection_f64_0_2.from_vec(a_i)
1533 elif order == 4:
1534 a_k = shamrock.math.SymTensorCollection_f64_0_3.from_vec(a_i)
1535 elif order == 5:
1536 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_i)
1537 else:
1538 raise ValueError("Invalid order")
1539
1540 if do_print:
1541 print("a_k =", a_k)
1542
1543 if order == 1:
1544 a_kp = shamrock.math.SymTensorCollection_f64_0_0.from_vec(a_ip)
1545 elif order == 2:
1546 a_kp = shamrock.math.SymTensorCollection_f64_0_1.from_vec(a_ip)
1547 elif order == 3:
1548 a_kp = shamrock.math.SymTensorCollection_f64_0_2.from_vec(a_ip)
1549 elif order == 4:
1550 a_kp = shamrock.math.SymTensorCollection_f64_0_3.from_vec(a_ip)
1551 elif order == 5:
1552 a_kp = shamrock.math.SymTensorCollection_f64_0_4.from_vec(a_ip)
1553 else:
1554 raise ValueError("Invalid order")
1555
1556 if do_print:
1557 print("a_kp =", a_kp)
1558
1559 if order == 1:
1560 result = shamrock.phys.contract_grav_moment_to_force_1(a_k, dM_k)
1561 elif order == 2:
1562 result = shamrock.phys.contract_grav_moment_to_force_2(a_k, dM_k)
1563 elif order == 3:
1564 result = shamrock.phys.contract_grav_moment_to_force_3(a_k, dM_k)
1565 elif order == 4:
1566 result = shamrock.phys.contract_grav_moment_to_force_4(a_k, dM_k)
1567 elif order == 5:
1568 result = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
1569 else:
1570 raise ValueError("Invalid order")
1571
1572 Gconst = 1 # let's just set the grav constant to 1
1573 force_i = -Gconst * np.array(result)
1574 if do_print:
1575 print("force_i =", force_i)
1576
1577 if order == 1:
1578 result_offset = shamrock.phys.contract_grav_moment_to_force_1(a_kp, dM_k_offset)
1579 elif order == 2:
1580 result_offset = shamrock.phys.contract_grav_moment_to_force_2(a_kp, dM_k_offset)
1581 elif order == 3:
1582 result_offset = shamrock.phys.contract_grav_moment_to_force_3(a_kp, dM_k_offset)
1583 elif order == 4:
1584 result_offset = shamrock.phys.contract_grav_moment_to_force_4(a_kp, dM_k_offset)
1585 elif order == 5:
1586 result_offset = shamrock.phys.contract_grav_moment_to_force_5(a_kp, dM_k_offset)
1587 else:
1588 raise ValueError("Invalid order")
1589
1590 force_i_offset = -Gconst * np.array(result_offset)
1591 if do_print:
1592 print("force_i_offset =", force_i_offset)
1593
1594 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
1595 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1596 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
1597
1598 angle = (b_A_size + b_B_size) / b_dist
1599
1600 delta_A = np.linalg.norm(np.array(s_A) - np.array(s_Ap))
1601
1602 if do_print:
1603 print("b_A_size =", b_A_size)
1604 print("b_B_size =", b_B_size)
1605 print("b_dist =", b_dist)
1606 print("angle =", angle)
1607
1608 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 ;) ).
1615 plt.figure()
1616 for order in range(1, 6):
1617 print("--------------------------------")
1618 print(f"Running FMM order = {order}")
1619 print("--------------------------------")
1620
1621 # set seed
1622 rng = np.random.default_rng(111)
1623
1624 N = 50000
1625
1626 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
1627 x_i_all = []
1628 for i in range(N):
1629 x_i_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1630
1631 # same for x_j
1632 x_j_all = []
1633 for i in range(N):
1634 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1635
1636 box_scale_fact_all = np.linspace(0, 0.1, N).tolist()
1637
1638 # same for box_1_center
1639 s_A_all = []
1640 s_Ap_all = []
1641 for p, box_scale_fact in zip(x_i_all, box_scale_fact_all):
1642 s_A_all.append(
1643 (
1644 p[0] + box_scale_fact * rng.uniform(-1, 1),
1645 p[1] + box_scale_fact * rng.uniform(-1, 1),
1646 p[2] + box_scale_fact * rng.uniform(-1, 1),
1647 )
1648 )
1649 s_Ap_all.append(
1650 (
1651 p[0] + box_scale_fact * rng.uniform(-1, 1),
1652 p[1] + box_scale_fact * rng.uniform(-1, 1),
1653 p[2] + box_scale_fact * rng.uniform(-1, 1),
1654 )
1655 )
1656
1657 # same for box_2_center
1658 s_B_all = []
1659 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1660 s_B_all.append(
1661 (
1662 p[0] + box_scale_fact * rng.uniform(-1, 1),
1663 p[1] + box_scale_fact * rng.uniform(-1, 1),
1664 p[2] + box_scale_fact * rng.uniform(-1, 1),
1665 )
1666 )
1667
1668 angles = []
1669 delta_A_all = []
1670 rel_errors = []
1671
1672 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):
1673 force_i, force_i_offset, angle, delta_A = test_grav_moment_offset(
1674 x_i, x_j, s_A, s_Ap, s_B, m_j, order, do_print=False
1675 )
1676
1677 abs_error = np.linalg.norm(force_i_offset - force_i)
1678 rel_error = abs_error / np.linalg.norm(force_i)
1679
1680 b_A_size = np.linalg.norm(np.array(s_A) - np.array(x_i))
1681 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1682 b_dist = np.linalg.norm(np.array(s_A) - np.array(s_B))
1683 angle = (b_A_size + b_B_size) / b_dist
1684
1685 if angle > 5.0 or angle < 1e-4:
1686 continue
1687
1688 angles.append(angle)
1689 delta_A_all.append(delta_A)
1690 rel_errors.append(rel_error)
1691
1692 print(f"Computed for {len(angles)} cases")
1693
1694 plt.scatter(angles, rel_errors, s=1, label=f"FMM order = {order}")
1695
1696
1697 plt.xlabel("Angle")
1698 plt.ylabel("$|f_{\\rm fmm} - f_{\\rm fmm, offset}| / |f_{\\rm fmm}|$ (Relative error)")
1699 plt.xscale("log")
1700 plt.yscale("log")
1701 plt.title("Grav moment translation error")
1702 plt.legend(loc="lower right")
1703 plt.grid()
1704 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)
)
1770 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1771 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)
1781 dM_k = shamrock.phys.get_dM_mat_5(D_n, Q_n_B)
1782 print("dM_k =", dM_k)
1783 a_k = shamrock.math.SymTensorCollection_f64_0_4.from_vec((0, 0, 0))
1784 print("a_k =", a_k)
1785 force_i_fmm_sA_null = shamrock.phys.contract_grav_moment_to_force_5(a_k, dM_k)
1786 Gconst = 1 # let's just set the grav constant to 1
1787 force_i_fmm_sA_null = -Gconst * np.array(force_i_fmm_sA_null)
1788 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
1793 def analytic_force_i(x_i, x_j, Gconst):
1794 force_i_direct = (x_j[0] - x_i[0], x_j[1] - x_i[1], x_j[2] - x_i[2])
1795 force_i_direct /= np.linalg.norm(force_i_direct) ** 3
1796 force_i_direct *= m_i
1797 return force_i_direct
1798
1799
1800 force_i_exact = analytic_force_i(x_i, x_j, Gconst)
1801 print("force_i_exact =", force_i_exact)
force_i_exact = [-1. 0. 0.]
1804 abs_error = np.linalg.norm(force_i - force_i_exact)
1805 rel_error = abs_error / np.linalg.norm(force_i)
1806
1807
1808 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1809 b_dist = np.linalg.norm(np.array(x_i) - np.array(s_B))
1810 angle = (b_B_size) / b_dist
1811
1812 print("abs_error =", abs_error)
1813 print("rel_error =", rel_error)
1814 print("angle =", angle)
1815
1816 assert rel_error < 1e-2
abs_error = 0.0014495314137263013
rel_error = 0.001447602143490707
angle = 0.23249527748763862
Let’s code MM in a box
1822 def run_mm(x_i, x_j, s_B, m_j, order, do_print):
1823 b_j = (x_j[0] - s_B[0], x_j[1] - s_B[1], x_j[2] - s_B[2])
1824 r = (s_B[0] - x_i[0], s_B[1] - x_i[1], s_B[2] - x_i[2])
1825
1826 if do_print:
1827 print("x_i =", x_i)
1828 print("x_j =", x_j)
1829 print("s_B =", s_B)
1830 print("b_j =", b_j)
1831 print("r =", r)
1832
1833 # compute the tensor product of the displacment
1834 if order == 1:
1835 Q_n_B = shamrock.math.SymTensorCollection_f64_0_0.from_vec(b_j)
1836 elif order == 2:
1837 Q_n_B = shamrock.math.SymTensorCollection_f64_0_1.from_vec(b_j)
1838 elif order == 3:
1839 Q_n_B = shamrock.math.SymTensorCollection_f64_0_2.from_vec(b_j)
1840 elif order == 4:
1841 Q_n_B = shamrock.math.SymTensorCollection_f64_0_3.from_vec(b_j)
1842 elif order == 5:
1843 Q_n_B = shamrock.math.SymTensorCollection_f64_0_4.from_vec(b_j)
1844 else:
1845 raise ValueError("Invalid order")
1846
1847 # multiply by mass to get the mass moment
1848 Q_n_B *= m_j
1849
1850 if do_print:
1851 print("Q_n_B =", Q_n_B)
1852
1853 # green function gradients
1854 if order == 1:
1855 D_n = shamrock.phys.green_func_grav_cartesian_1_1(r)
1856 elif order == 2:
1857 D_n = shamrock.phys.green_func_grav_cartesian_1_2(r)
1858 elif order == 3:
1859 D_n = shamrock.phys.green_func_grav_cartesian_1_3(r)
1860 elif order == 4:
1861 D_n = shamrock.phys.green_func_grav_cartesian_1_4(r)
1862 elif order == 5:
1863 D_n = shamrock.phys.green_func_grav_cartesian_1_5(r)
1864 else:
1865 raise ValueError("Invalid order")
1866
1867 if do_print:
1868 print("D_n =", D_n)
1869
1870 if order == 1:
1871 result = shamrock.phys.contract_grav_moment_to_force_1(Q_n_B, D_n)
1872 elif order == 2:
1873 result = shamrock.phys.contract_grav_moment_to_force_2(Q_n_B, D_n)
1874 elif order == 3:
1875 result = shamrock.phys.contract_grav_moment_to_force_3(Q_n_B, D_n)
1876 elif order == 4:
1877 result = shamrock.phys.contract_grav_moment_to_force_4(Q_n_B, D_n)
1878 elif order == 5:
1879 result = shamrock.phys.contract_grav_moment_to_force_5(Q_n_B, D_n)
1880 else:
1881 raise ValueError("Invalid order")
1882
1883 Gconst = 1 # let's just set the grav constant to 1
1884 force_i = -Gconst * np.array(result)
1885 if do_print:
1886 print("force_i =", force_i)
1887
1888 force_i_exact = analytic_force_i(x_i, x_j, Gconst)
1889 if do_print:
1890 print("force_i_exact =", force_i_exact)
1891
1892 b_B_size = np.linalg.norm(np.array(s_B) - np.array(x_j))
1893 b_dist = np.linalg.norm(np.array(x_i) - np.array(s_B))
1894
1895 angle = (b_B_size) / b_dist
1896
1897 if do_print:
1898 print("b_A_size =", b_A_size)
1899 print("b_B_size =", b_B_size)
1900 print("b_dist =", b_dist)
1901 print("angle =", angle)
1902
1903 return force_i, force_i_exact, angle
1908 plt.figure()
1909 for order in range(1, 6):
1910 print("--------------------------------")
1911 print(f"Running MM order = {order}")
1912 print("--------------------------------")
1913
1914 # set seed
1915 rng = np.random.default_rng(111)
1916
1917 N = 50000
1918
1919 # generate a random set of position in a box of bounds (-1,1)x(-1,1)x(-1,1)
1920 x_i_all = []
1921 for i in range(N):
1922 x_i_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1923
1924 # same for x_j
1925 x_j_all = []
1926 for i in range(N):
1927 x_j_all.append((rng.uniform(-1, 1), rng.uniform(-1, 1), rng.uniform(-1, 1)))
1928
1929 box_scale_fact_all = np.linspace(0, 0.1, N).tolist()
1930
1931 # same for box_2_center
1932 s_B_all = []
1933 for p, box_scale_fact in zip(x_j_all, box_scale_fact_all):
1934 s_B_all.append(
1935 (
1936 p[0] + box_scale_fact * rng.uniform(-1, 1),
1937 p[1] + box_scale_fact * rng.uniform(-1, 1),
1938 p[2] + box_scale_fact * rng.uniform(-1, 1),
1939 )
1940 )
1941
1942 angles = []
1943 rel_errors = []
1944
1945 for x_i, x_j, s_B in zip(x_i_all, x_j_all, s_B_all):
1946 force_i, force_i_exact, angle = run_mm(x_i, x_j, s_B, m_j, order, do_print=False)
1947
1948 abs_error = np.linalg.norm(force_i - force_i_exact)
1949 rel_error = abs_error / np.linalg.norm(force_i_exact)
1950
1951 if angle > 5.0 or angle < 1e-4:
1952 continue
1953
1954 angles.append(angle)
1955 rel_errors.append(rel_error)
1956
1957 print(f"Computed for {len(angles)} cases")
1958
1959 plt.scatter(angles, rel_errors, s=1, label=f"MM order = {order}")
1960
1961
1962 def plot_powerlaw(order, center_y):
1963 X = [1e-3, 1e-2 / 3, 1e-1]
1964 Y = [center_y * (x) ** order for x in X]
1965 plt.plot(X, Y, linestyle="dashed", color="black")
1966 bbox = dict(boxstyle="round", fc="blanchedalmond", ec="orange", alpha=0.9)
1967 plt.text(X[1], Y[1], f"$\\propto x^{order}$", fontsize=9, bbox=bbox)
1968
1969
1970 plot_powerlaw(1, 1)
1971 plot_powerlaw(2, 1)
1972 plot_powerlaw(3, 1)
1973 plot_powerlaw(4, 1)
1974 plot_powerlaw(5, 1)
1975
1976 plt.xlabel("Angle")
1977 plt.ylabel("Relative Error")
1978 plt.xscale("log")
1979 plt.yscale("log")
1980 plt.title("MM precision")
1981 plt.legend(loc="lower right")
1982 plt.grid()
1983 plt.show()
--------------------------------
Running MM order = 1
--------------------------------
Computed for 49923 cases
--------------------------------
Running MM order = 2
--------------------------------
Computed for 49923 cases
--------------------------------
Running MM order = 3
--------------------------------
Computed for 49923 cases
--------------------------------
Running MM order = 4
--------------------------------
Computed for 49923 cases
--------------------------------
Running MM order = 5
--------------------------------
Computed for 49923 cases
Total running time of the script: (0 minutes 52.386 seconds)
Estimated memory usage: 193 MB