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