Select Git revision
multipole_accept.h
Matthieu Schaller authored
multipole_accept.h 12.32 KiB
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_MULTIPOLE_ACCEPT_H
#define SWIFT_MULTIPOLE_ACCEPT_H
/* Config parameters. */
#include <config.h>
/* Local includes */
#include "binomial.h"
#include "gravity_properties.h"
#include "integer_power.h"
#include "kernel_long_gravity.h"
#include "minmax.h"
#include "multipole_struct.h"
/**
* @brief Compute the inverse of the force estimator entering the MAC
*
* Note that in the unsofted case, the first condition is naturally
* never reached (as H == 0). In the non-periodic (non-truncated) case
* the second condition is never reached (as r_s == inf, r_s_inv == 0).
*
* @param H The spline softening length.
* @param r_s_inv The inverse of the scale of the gravity mesh.
* @param r2 The square of the distance between the multipoles.
*/
__attribute__((const)) INLINE static float gravity_f_MAC_inverse(
const float H, const float r_s_inv, const float r2) {
if (r2 < (25.f / 81.f) * H * H) {
/* Below softening radius */
return (25.f / 81.f) * H * H;
} else if (r_s_inv * r_s_inv * r2 > (25.f / 9.f)) {
/* Above truncation radius */
return (9.f / 25.f) * r_s_inv * r_s_inv * r2 * r2;
} else {
/* Normal Newtonian case */
return r2;
}
}
/**
* @brief Checks whether The multipole in B can be used to update the field
* tensor in A.
*
* We use the MAC of Dehnen 2014 eq. 16.
*
* Note: this is *not* symmetric in A<->B unless the purely geometric criterion * is used.
*
* @param props The properties of the gravity scheme.
* @param A The gravity tensors that we want to update (sink).
* @param B The gravity tensors that act as a source.
* @param r2 The square of the distance between the centres of mass of A and B.
* @param use_rebuild_sizes Are we considering the sizes at the last tree-build
* (1) or current sizes (0)?
* @param periodic Are we using periodic BCs?
*/
__attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
const struct gravity_props *props, const struct gravity_tensors *restrict A,
const struct gravity_tensors *restrict B, const float r2,
const int use_rebuild_sizes, const int periodic) {
/* Order of the expansion */
const int p = 2;
/* Sizes of the multipoles */
const float rho_A = use_rebuild_sizes ? A->r_max_rebuild : A->r_max;
const float rho_B = use_rebuild_sizes ? B->r_max_rebuild : B->r_max;
/* Max size of both multipoles */
const float rho_max = max(rho_A, rho_B);
/* Get the softening */
const float max_softening =
max(A->m_pole.max_softening, B->m_pole.max_softening);
/* Compute the error estimator (without the 1/M_B term that cancels out) */
float E_BA_term = 0.f;
for (int n = 0; n <= p; ++n) {
E_BA_term +=
binomial(p, n) * B->m_pole.power[n] * integer_powf(rho_A, p - n);
}
E_BA_term *= 8.f;
if (rho_A + rho_B > 0.f) {
E_BA_term *= rho_max;
;
E_BA_term /= (rho_A + rho_B);
}
/* Compute r^p = (r^2)^(p/2) */
const float r_to_p = integer_powf(r2, (p / 2));
float f_MAC_inv;
if (periodic && props->consider_truncation_in_MAC) {
f_MAC_inv = gravity_f_MAC_inverse(max_softening, props->r_s_inv, r2);
} else {
f_MAC_inv = r2;
}
/* Get the mimimal acceleration in A */
const float min_a_grav = A->m_pole.min_old_a_grav_norm;
/* Maximal mass */
const float M_max = max(A->m_pole.M_000, B->m_pole.M_000);
/* Get the relative tolerance */
const float eps = props->adaptive_tolerance;
/* Get the basic geometric critical angle */
const float theta_crit = props->theta_crit;
const float theta_crit2 = theta_crit * theta_crit;
/* Get the sum of the multipole sizes */
const float rho_sum = rho_A + rho_B;
if (props->use_advanced_MAC && props->use_gadget_tolerance) {
/* Gadget 4 paper -- eq. 36 */
const int power = SELF_GRAVITY_MULTIPOLE_ORDER - 1;
const float ratio = integer_powf(rho_max / sqrtf(r2), power);
const int cond_1 = M_max * ratio < eps * min_a_grav * f_MAC_inv;
const int cond_2 =
props->use_tree_below_softening || max_softening * max_softening < r2;
return cond_1 && cond_2;
} else if (props->use_advanced_MAC && !props->use_gadget_tolerance) {
#ifdef SWIFT_DEBUG_CHECKS
if (min_a_grav == 0.) error("Acceleration is 0");
#endif
/* Test the different conditions */
/* Condition 1: We are in the converging part of the Taylor expansion */
const int cond_1 = rho_sum * rho_sum < r2;
/* Condition 2: We are not below softening */
const int cond_2 =
props->use_tree_below_softening || max_softening * max_softening < r2;
/* Condition 3: The contribution is accurate enough
* (E_BA * (1 / r^(p)) * ((1 / r^2) * W) < eps * a_min) */
const int cond_3 = E_BA_term < eps * min_a_grav * r_to_p * f_MAC_inv;
return cond_1 && cond_2 && cond_3;
} else {
/* Condition 1: We are obeying the purely geometric criterion */
const int cond_1 = rho_sum * rho_sum < theta_crit2 * r2;
/* Condition 2: We are not below softening */
const int cond_2 =
props->use_tree_below_softening || max_softening * max_softening < r2;
return cond_1 && cond_2;
}
}
/**
* @brief Checks whether The multipole in B can be used to update the field
* tensor in A and whether the multipole in A can be used to update the field
* tensor in B.
*
* We use the MAC of Dehnen 2014 eq. 16.
*
* @param props The properties of the gravity scheme.
* @param A The first set of multipole and gravity tensors.
* @param B The second set of multipole and gravity tensors.
* @param r2 The square of the distance between the centres of mass of A and B.
* @param use_rebuild_sizes Are we considering the sizes at the last tree-build
* (1) or current sizes (0)?
* @param periodic Are we using periodic BCs?
*/
__attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept_symmetric(
const struct gravity_props *props, const struct gravity_tensors *restrict A,
const struct gravity_tensors *restrict B, const float r2,
const int use_rebuild_sizes, const int periodic) {
return gravity_M2L_accept(props, A, B, r2, use_rebuild_sizes, periodic) &&
gravity_M2L_accept(props, B, A, r2, use_rebuild_sizes, periodic);
}
/**
* Compute the distance above which an M2L kernel is allowed to be used. *
* This uses conservative assumptions to guarantee that all the possible cell
* pair interactions that need a direct interaction are below this distance.
*
* @param props The properties of the gravity scheme.
* @param size The size of the multipoles (here the cell size).
* @param max_softening The maximal softening accross all particles.
* @param min_a_grav The minimal acceleration accross all particles.
* @param max_mpole_power The maximum multipole power accross all the
* multipoles.
* @param periodic Are we using periodic BCs?
*/
__attribute__((nonnull, pure)) INLINE static float
gravity_M2L_min_accept_distance(
const struct gravity_props *props, const float size,
const float max_softening, const float min_a_grav,
const float max_mpole_power[SELF_GRAVITY_MULTIPOLE_ORDER + 1],
const int periodic) {
/* Order of the expansion */
const int p = SELF_GRAVITY_MULTIPOLE_ORDER;
float E_BA_term = 0.f;
for (int n = 0; n <= p; ++n) {
E_BA_term +=
binomial(p, n) * max_mpole_power[n] * integer_powf(size, p - n);
}
E_BA_term *= 4.f;
/* Get the basic geometric critical angle */
const float theta_crit = props->theta_crit;
const float theta_crit2 = theta_crit * theta_crit;
/* Get the sum of the multipole sizes */
const float size_sum = 2. * size;
/* Get the relative tolerance */
const float eps = props->adaptive_tolerance;
if (props->use_advanced_MAC) {
/* Distance obtained by solving for the geometric criterion with theta = 1
*/
const float dist_tree = size_sum;
const float dist_adapt =
powf(E_BA_term / (eps * min_a_grav), 1.f / (p + 2.f));
/* Distance obtained by demanding > softening */
const float dist_soft =
props->use_tree_below_softening ? 0.f : max_softening;
return max3(dist_tree, dist_adapt, dist_soft);
} else {
/* Distance obtained by solving for the geometric criterion */
const float dist_tree = sqrtf(size_sum * size_sum / theta_crit2);
/* Distance obtained by demanding > softening */
const float dist_soft =
props->use_tree_below_softening ? 0.f : max_softening;
return max(dist_tree, dist_soft);
}
}
/**
* @brief Checks whether The multipole in B can be used to update the particle
* pa *
* We use the MAC of Dehnen 2014 eq. 16.
*
* @param props The properties of the gravity scheme.
* @param pa The particle we want to compute forces for (sink)
* @param B The gravity tensors that act as a source.
* @param r2 The square of the distance between pa and the centres of mass of B.
* @param periodic Are we using periodic BCs?
*/
__attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
const struct gravity_props *props, const struct gpart *pa,
const struct gravity_tensors *B, const float r2, const int periodic) {
/* Order of the expansion */
const int p = 2;
/* Sizes of the multipoles */
const float rho_B = B->r_max;
/* Get the maximal softening */
const float max_softening =
max(B->m_pole.max_softening, gravity_get_softening(pa, props));
#ifdef SWIFT_DEBUG_CHECKS
if (rho_B == 0.) error("Size of multipole B is 0!");
#endif
/* Compute the error estimator (without the 1/M_B term that cancels out) */
const float E_BA_term = 8.f * B->m_pole.power[p];
/* Compute r^p = (r^2)^(p/2) */
const float r_to_p = integer_powf(r2, (p / 2));
float f_MAC_inv;
if (periodic && props->consider_truncation_in_MAC) {
f_MAC_inv = gravity_f_MAC_inverse(max_softening, props->r_s_inv, r2);
} else {
f_MAC_inv = r2;
}
/* Get the estimate of the acceleration */
const float old_a_grav = pa->old_a_grav_norm;
/* Get the relative tolerance */
const float eps = props->adaptive_tolerance;
/* Get the basic geometric critical angle */
const float theta_crit = props->theta_crit;
const float theta_crit2 = theta_crit * theta_crit;
if (props->use_advanced_MAC && props->use_gadget_tolerance) {
/* Gadget 4 paper -- eq. 12 */
const int power = SELF_GRAVITY_MULTIPOLE_ORDER;
const float ratio = integer_powf(rho_B / sqrtf(r2), power);
const int cond_1 = B->m_pole.M_000 * ratio < eps * old_a_grav * f_MAC_inv;
const int cond_2 =
props->use_tree_below_softening || max_softening * max_softening < r2;
return cond_1 && cond_2;
} else if (props->use_advanced_MAC && !props->use_gadget_tolerance) {
#ifdef SWIFT_DEBUG_CHECKS
if (old_a_grav == 0.) error("Acceleration is 0");
#endif
/* Test the different conditions */
/* Condition 1: We are in the converging part of the Taylor expansion */
const int cond_1 = rho_B * rho_B < r2;
/* Condition 2: We are not below softening */
const int cond_2 =
props->use_tree_below_softening || max_softening * max_softening < r2;
/* Condition 3: The contribution is accurate enough
* (E_BA * (1 / r^(p)) * ((1 / r^2) * W) < eps * a_min) */
const int cond_3 = E_BA_term < eps * old_a_grav * r_to_p * f_MAC_inv;
return cond_1 && cond_2 && cond_3;
} else {
/* Condition 1: We are obeying the purely geometric criterion */
const int cond_1 = rho_B * rho_B < theta_crit2 * r2;
/* Condition 2: We are not below softening */
const int cond_2 =
props->use_tree_below_softening || max_softening * max_softening < r2;
return cond_1 && cond_2;
}
}
#endif /* SWIFT_MULTIPOLE_ACCEPT_H */