Skip to content
Snippets Groups Projects
Commit a651fa11 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Implemented the new updated MAC

parent 69f724ab
No related branches found
No related tags found
1 merge request!1077Improved multipole acceptance criterion (MAC)
...@@ -32,6 +32,7 @@ ...@@ -32,6 +32,7 @@
#include "gravity.h" #include "gravity.h"
#include "kernel_gravity.h" #include "kernel_gravity.h"
#include "kernel_long_gravity.h" #include "kernel_long_gravity.h"
#include "restart.h"
#define gravity_props_default_a_smooth 1.25f #define gravity_props_default_a_smooth 1.25f
#define gravity_props_default_r_cut_max 4.5f #define gravity_props_default_r_cut_max 4.5f
...@@ -86,8 +87,6 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params, ...@@ -86,8 +87,6 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
/* Opening angle */ /* Opening angle */
p->theta_crit = parser_get_param_double(params, "Gravity:theta"); p->theta_crit = parser_get_param_double(params, "Gravity:theta");
if (p->theta_crit >= 1.) error("Theta too large. FMM won't converge."); if (p->theta_crit >= 1.) error("Theta too large. FMM won't converge.");
p->theta_crit2 = p->theta_crit * p->theta_crit;
p->theta_crit_inv = 1. / p->theta_crit;
/* Mesh dithering */ /* Mesh dithering */
if (periodic && !with_external_potential) { if (periodic && !with_external_potential) {
......
...@@ -26,10 +26,6 @@ ...@@ -26,10 +26,6 @@
#include <hdf5.h> #include <hdf5.h>
#endif #endif
/* Local includes. */
#include "kernel_gravity.h"
#include "restart.h"
/* Forward declarations */ /* Forward declarations */
struct cosmology; struct cosmology;
struct phys_const; struct phys_const;
...@@ -58,14 +54,14 @@ struct gravity_props { ...@@ -58,14 +54,14 @@ struct gravity_props {
/* -------------- Properties of the FFM gravity ---------------------- */ /* -------------- Properties of the FFM gravity ---------------------- */
/*! Tree opening angle (Multipole acceptance criterion) */ /*! Are we using the adaptive opening angle? */
double theta_crit; int use_adaptive_tolerance;
/*! Square of opening angle */ /*! Accuracy parameter of the advanced MAC */
double theta_crit2; float adaptive_tolerance;
/*! Inverse of opening angle */ /*! Tree opening angle (Multipole acceptance criterion) */
double theta_crit_inv; double theta_crit;
/* ------------- Properties of the softened gravity ------------------ */ /* ------------- Properties of the softened gravity ------------------ */
......
...@@ -23,67 +23,160 @@ ...@@ -23,67 +23,160 @@
#include "../config.h" #include "../config.h"
/* Local includes */ /* Local includes */
#include "binomial.h"
#include "integer_power.h"
#include "minmax.h"
#include "multipole_struct.h" #include "multipole_struct.h"
/** /**
* @brief Checks whether a cell-cell interaction can be appromixated by a M-M * @brief Checks whether The multipole in B can be used to update the field
* interaction using the distance and cell radius. * tensor in A.
* *
* We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179, * We use the MAC of Dehnen 2014 eq. 16.
* Issue 1, pp.27-42, equation 10.
* *
* We also additionally check that the distance between the multipoles * Note: this is *not* symmetric in A<->B unless the purely geometric criterion
* is larger than the softening lengths (here the distance at which * is used.
* the gravity becomes Newtonian again, not the Plummer-equivalent quantity).
* *
* @param r_crit_a The size of the multipole A. * @param props The properties of the gravity scheme.
* @param r_crit_b The size of the multipole B. * @param A The gravity tensors that we want to update (sink).
* @param theta_crit2 The square of the critical opening angle. * @param B The gravity tensors that act as a source.
* @param r2 Square of the distance (periodically wrapped) between the * @param r2 The square of the distance between the centres of mass of A and B.
* multipoles.
* @param epsilon_a The maximal softening length of any particle in A.
* @param epsilon_b The maximal softening length of any particle in B.
*/ */
__attribute__((always_inline, const)) INLINE static int gravity_M2L_accept( __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
const double r_crit_a, const double r_crit_b, const double theta_crit2, const struct gravity_props *props, const struct gravity_tensors *restrict A,
const double r2, const double epsilon_a, const double epsilon_b) { const struct gravity_tensors *restrict B, const float r2) {
const double size = r_crit_a + r_crit_b; /* Order of the expansion */
const double size2 = size * size; const int p = SELF_GRAVITY_MULTIPOLE_ORDER;
const double epsilon_a2 = epsilon_a * epsilon_a;
const double epsilon_b2 = epsilon_b * epsilon_b;
// MATTHIEU: Make this mass-dependent ? /* 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(A->r_max, p - n);
}
E_BA_term *= 8.f;
E_BA_term *= max(A->r_max, B->r_max);
E_BA_term /= (A->r_max + B->r_max);
/* Multipole acceptance criterion (Dehnen 2002, eq.10) */ /* Compute r^(p+2) */
return (r2 * theta_crit2 > size2) && (r2 > epsilon_a2) && (r2 > epsilon_b2); #if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
#else
const float r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
#endif
/* Get the mimimal acceleration in A */
const float min_a_grav = A->m_pole.min_old_a_grav_norm;
/* Get the maximal softening length in B */
const float max_softening = B->m_pole.max_softening;
/* 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 = A->r_max + B->r_max;
if (props->use_adaptive_tolerance) {
/* 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 = max_softening * max_softening < r2;
/* Condition 3: The contribution is accurate enough
* (E_BA / r^(p+2) < eps a_min) */
const int cond_3 = E_BA_term < eps * min_a_grav * r_to_p_plus2;
return cond_1 && cond_2 && cond_3;
} else {
/* Condition 1: We are in the converging part of the Taylor expansion */
const int cond_1 = rho_sum * rho_sum < theta_crit2 * r2;
/* Condition 2: We are not below softening */
const int cond_2 = max_softening * max_softening < r2;
return cond_1 && cond_2;
}
} }
/** /**
* @brief Checks whether a particle-cell interaction can be appromixated by a * @brief Checks whether The multipole in B can be used to update the particle
* M2P interaction using the distance and cell radius. * pa
*
* We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
* Issue 1, pp.27-42, equation 10.
* *
* We also additionally check that the distance between the particle and the * We use the MAC of Dehnen 2014 eq. 16.
* multipole is larger than the softening length (here the distance at which
* the gravity becomes Newtonian again, not the Plummer-equivalent quantity).
* *
* @param r_max2 The square of the size of the multipole. * @param props The properties of the gravity scheme.
* @param theta_crit2 The square of the critical opening angle. * @param pa The particle we want to compute forces for (sink)
* @param r2 Square of the distance (periodically wrapped) between the * @param B The gravity tensors that act as a source.
* particle and the multipole. * @param r2 The square of the distance between pa and the centres of mass of B.
* @param epsilon The softening length of the particle.
*/ */
__attribute__((always_inline, const)) INLINE static int gravity_M2P_accept( __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
const float r_max2, const float theta_crit2, const float r2, const struct gravity_props *props, const struct gpart *pa,
const float epsilon) { const struct gravity_tensors *B, const float r2) {
/* Order of the expansion */
const int p = SELF_GRAVITY_MULTIPOLE_ORDER;
/* Compute the error estimator (without the 1/M_B term that cancels out) */
float E_BA_term = 8.f * B->m_pole.power[p];
/* Compute r^(p+2) */
#if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
#else
const float r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
#endif
/* Get the estimate of the acceleration */
const float old_a_grav = pa->old_a_grav_norm;
/* Get the maximal softening length in B */
const float max_softening = B->m_pole.max_softening;
/* 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_adaptive_tolerance) {
/* Test the different conditions */
/* Condition 1: We are in the converging part of the Taylor expansion */
const int cond_1 = (B->r_max) * (B->r_max) < r2;
/* Condition 2: We are not below softening */
const int cond_2 = max_softening * max_softening < r2;
/* Condition 3: The contribution is accurate enough
* (E_BA / r^(p+2) < eps * a) */
const int cond_3 = E_BA_term < eps * old_a_grav * r_to_p_plus2;
return cond_1 && cond_2 && cond_3;
} else {
/* Condition 1: We are in the converging part of the Taylor expansion */
const int cond_1 = (B->r_max) * (B->r_max) < theta_crit2 * r2;
// MATTHIEU: Make this mass-dependent ? /* Condition 2: We are not below softening */
const int cond_2 = max_softening * max_softening < r2;
/* Multipole acceptance criterion (Dehnen 2002, eq.10) */ return cond_1 && cond_2;
return (r2 * theta_crit2 > r_max2) && (r2 > epsilon * epsilon); }
} }
#endif /* SWIFT_MULTIPOLE_ACCEPT_H */ #endif /* SWIFT_MULTIPOLE_ACCEPT_H */
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment