Commit 8446277c authored by Matthieu Schaller's avatar Matthieu Schaller

Also clarified and simplied the expressions for D_tilde in the long-range truncated case

parent 6f6c3499
......@@ -303,35 +303,55 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
Dt_1 = derivs.chi_0 * r_inv;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
const float r_inv2 = r_inv * r_inv;
Dt_2 = (r * derivs.chi_1 - derivs.chi_0) * r_inv2;
/* -chi^0 r_i^2 + chi^1 r_i^1 */
Dt_2 = derivs.chi_1 - derivs.chi_0 * r_inv;
Dt_2 = Dt_2 * r_inv;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
const float r_inv3 = r_inv2 * r_inv;
Dt_3 =
(r * r * derivs.chi_2 - 3.f * r * derivs.chi_1 + 3.f * derivs.chi_0) *
r_inv3;
/* 3chi^0 r_i^3 - 3 chi^1 r_i^2 + chi^2 r_i^1 */
Dt_3 = derivs.chi_0 * r_inv - derivs.chi_1;
Dt_3 = Dt_3 * 3.f;
Dt_3 = Dt_3 * r_inv + derivs.chi_2;
Dt_3 = Dt_3 * r_inv;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
const float r_inv4 = r_inv3 * r_inv;
Dt_4 = (r * r * r * derivs.chi_3 - 6.f * r * r * derivs.chi_2 +
15.f * r * derivs.chi_1 - 15.f * derivs.chi_0) *
r_inv4;
/* -15chi^0 r_i^4 + 15 chi^1 r_i^3 - 6 chi^2 r_i^2 + chi^3 r_i^1 */
Dt_4 = -derivs.chi_0 * r_inv + derivs.chi_1;
Dt_4 = Dt_4 * 15.f;
Dt_4 = Dt_4 * r_inv - 6.f * derivs.chi_2;
Dt_4 = Dt_4 * r_inv + derivs.chi_3;
Dt_4 = Dt_4 * r_inv;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
const float r_inv5 = r_inv4 * r_inv;
Dt_5 = (r * r * r * r * derivs.chi_4 - 10.f * r * r * r * derivs.chi_3 +
45.f * r * r * derivs.chi_2 - 105.f * r * derivs.chi_1 +
105.f * derivs.chi_0) *
r_inv5;
/* 105chi^0 r_i^5 - 105 chi^1 r_i^4 + 45 chi^2 r_i^3 - 10 chi^3 r_i^2 +
* chi^4 r_i^1 */
Dt_5 = derivs.chi_0 * r_inv - derivs.chi_1;
Dt_5 = Dt_5 * 105.f;
Dt_5 = Dt_5 * r_inv + 45.f * derivs.chi_2;
Dt_5 = Dt_5 * r_inv - 10.f * derivs.chi_3;
Dt_5 = Dt_5 * r_inv + derivs.chi_4;
Dt_5 = Dt_5 * r_inv;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
const float r_inv6 = r_inv5 * r_inv;
Dt_6 = (r * r * r * r * r * derivs.chi_5 -
15.f * r * r * r * r * derivs.chi_4 +
105.f * r * r * r * derivs.chi_3 - 420.f * r * r * derivs.chi_2 +
945.f * r * derivs.chi_1 - 945.f * derivs.chi_0) *
r_inv6;
/* -945chi^0 r_i^6 + 945 chi^1 r_i^5 - 420 chi^2 r_i^4 + 105 chi^3 r_i^3 -
* 15 chi^4 r_i^2 + chi^5 r_i^1 */
Dt_6 = -derivs.chi_0 * r_inv + derivs.chi_1;
Dt_6 = Dt_6 * 945.f;
Dt_6 = Dt_6 * r_inv - 420.f * derivs.chi_2;
Dt_6 = Dt_6 * r_inv + 105.f * derivs.chi_3;
Dt_6 = Dt_6 * r_inv - 15.f * derivs.chi_4;
Dt_6 = Dt_6 * r_inv + derivs.chi_5;
Dt_6 = Dt_6 * r_inv;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
......@@ -570,34 +590,53 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
Dt_1 = derivs.chi_0 * r_inv;
const float r_inv2 = r_inv * r_inv;
Dt_2 = (r * derivs.chi_1 - derivs.chi_0) * r_inv2;
/* -chi^0 r_i^2 + chi^1 r_i^1 */
Dt_2 = derivs.chi_1 - derivs.chi_0 * r_inv;
Dt_2 = Dt_2 * r_inv;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
const float r_inv3 = r_inv2 * r_inv;
Dt_3 =
(r * r * derivs.chi_2 - 3.f * r * derivs.chi_1 + 3.f * derivs.chi_0) *
r_inv3;
/* 3chi^0 r_i^3 - 3 chi^1 r_i^2 + chi^2 r_i^1 */
Dt_3 = derivs.chi_0 * r_inv - derivs.chi_1;
Dt_3 = Dt_3 * 3.f;
Dt_3 = Dt_3 * r_inv + derivs.chi_2;
Dt_3 = Dt_3 * r_inv;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
const float r_inv4 = r_inv3 * r_inv;
Dt_4 = (r * r * r * derivs.chi_3 - 6.f * r * r * derivs.chi_2 +
15.f * r * derivs.chi_1 - 15.f * derivs.chi_0) *
r_inv4;
/* -15chi^0 r_i^4 + 15 chi^1 r_i^3 - 6 chi^2 r_i^2 + chi^3 r_i^1 */
Dt_4 = -derivs.chi_0 * r_inv + derivs.chi_1;
Dt_4 = Dt_4 * 15.f;
Dt_4 = Dt_4 * r_inv - 6.f * derivs.chi_2;
Dt_4 = Dt_4 * r_inv + derivs.chi_3;
Dt_4 = Dt_4 * r_inv;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
const float r_inv5 = r_inv4 * r_inv;
Dt_5 = (r * r * r * r * derivs.chi_4 - 10.f * r * r * r * derivs.chi_3 +
45.f * r * r * derivs.chi_2 - 105.f * r * derivs.chi_1 +
105.f * derivs.chi_0) *
r_inv5;
/* 105chi^0 r_i^5 - 105 chi^1 r_i^4 + 45 chi^2 r_i^3 - 10 chi^3 r_i^2 +
* chi^4 r_i^1 */
Dt_5 = derivs.chi_0 * r_inv - derivs.chi_1;
Dt_5 = Dt_5 * 105.f;
Dt_5 = Dt_5 * r_inv + 45.f * derivs.chi_2;
Dt_5 = Dt_5 * r_inv - 10.f * derivs.chi_3;
Dt_5 = Dt_5 * r_inv + derivs.chi_4;
Dt_5 = Dt_5 * r_inv;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
const float r_inv6 = r_inv5 * r_inv;
Dt_6 = (r * r * r * r * r * derivs.chi_5 -
15.f * r * r * r * r * derivs.chi_4 +
105.f * r * r * r * derivs.chi_3 - 420.f * r * r * derivs.chi_2 +
945.f * r * derivs.chi_1 - 945.f * derivs.chi_0) *
r_inv6;
/* -945chi^0 r_i^6 + 945 chi^1 r_i^5 - 420 chi^2 r_i^4 + 105 chi^3 r_i^3 -
* 15 chi^4 r_i^2 + chi^5 r_i^1 */
Dt_6 = -derivs.chi_0 * r_inv + derivs.chi_1;
Dt_6 = Dt_6 * 945.f;
Dt_6 = Dt_6 * r_inv - 420.f * derivs.chi_2;
Dt_6 = Dt_6 * r_inv + 105.f * derivs.chi_3;
Dt_6 = Dt_6 * r_inv - 15.f * derivs.chi_4;
Dt_6 = Dt_6 * r_inv + derivs.chi_5;
Dt_6 = Dt_6 * r_inv;
#endif
}
......
......@@ -150,7 +150,21 @@ In the case $u<1$ and using $f(u)$ given by \ref{eq:fmm:potential}, we can simpl
\mathsf{\tilde{D}}_{5} &= (-315u^3 + 420u) \times H^{-5}, \nonumber \\
\mathsf{\tilde{D}}_{6} &= (315u^2 - 1260) \times H^{-6}. \nonumber
\end{align}
We can now write out all the derivatives used in the M2L and M2P kernels:
These expressions only use low powers of $u$ and, in particular, no terms
involving $1/u$ as would be the case when using a cubic spline kernel for
$f(u)$. This makes this choice of softening kernel much faster to evaluate than
ones using divisions. Similarly, the expressions in the periodic case for $u>1$
can be simplified to:
\begin{align}
\mathsf{\tilde{D}}_{1} &= \chi r^{-1}, \nonumber \\
\mathsf{\tilde{D}}_{2} &= -\chi r^{-2} + \chi' r^{-1}, \nonumber \\
\mathsf{\tilde{D}}_{3} &= 3\chi r^{-3} - 3\chi' r^{-2} + \chi'' r^{-1}, \nonumber \\
\mathsf{\tilde{D}}_{4} &= -15\chi r^{-4} + 15\chi' r^{-3} - 6\chi''r^{-2} + \chi^{(3)} r^{-1}, \nonumber \\
\mathsf{\tilde{D}}_{5} &= 105\chi r^{-5} -105\chi' r^{-4} + 45\chi''r^{-3} - 10\chi^{(3)} r^{-2} + \chi^{(4)} r^{-1}\nonumber, \\
\mathsf{\tilde{D}}_{6} &= -945\chi r^{-6} + 945 \chi' r^{-5} -420 \chi'' r^{-4} + 105 \chi^{(3)} r^{-3} - 15\chi^{(4)} r^{-2} + \chi^{(5)} r^{-1}. \nonumber
\end{align}
We can now write out all the derivatives used in the M2L and
M2P kernels:
\begin{align}
\mathsf{D}_{000}(\mathbf{r}) = \varphi (\mathbf{r}, r_s, H) =
\mathsf{\tilde{D}}_{1} \nonumber
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment