Skip to content

GitLab

Commit 23129ec0 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added short-range calculation of the potential.

parent 3b8cb588
Hide whitespace changes
Inline Side-by-side
src/gravity/Default/gravity_iact.h
View file @ 23129ec0
......@@ -40,7 +40,8 @@
* @param f_ij (return) The force intensity.
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij) {
float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij,
float *pot_ij) {
/* Get the inverse distance */
const float r_inv = 1.f / sqrtf(r2);
......@@ -50,17 +51,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
/* Get Newtonian gravity */
*f_ij = mass * r_inv * r_inv * r_inv;
*pot_ij = -mass * r_inv;
} else {
const float r = r2 * r_inv;
const float ui = r * h_inv;
float W_ij;
kernel_grav_eval(ui, &W_ij);
float W_f_ij, W_pot_ij;
kernel_grav_force_eval(ui, &W_f_ij);
kernel_grav_pot_eval(ui, &W_pot_ij);
/* Get softened gravity */
*f_ij = mass * h_inv3 * W_ij;
*f_ij = mass * h_inv3 * W_f_ij;
*pot_ij = mass * h_inv * W_pot_ij;
}
}
......@@ -81,7 +85,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
float r2, float h2, float h_inv, float h_inv3, float mass, float rlr_inv,
float *f_ij) {
float *f_ij, float *pot_ij) {
/* Get the inverse distance */
const float r_inv = 1.f / sqrtf(r2);
......@@ -92,23 +96,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
/* Get Newtonian gravity */
*f_ij = mass * r_inv * r_inv * r_inv;
*pot_ij = -mass * r_inv;
} else {
const float ui = r * h_inv;
float W_ij;
float W_f_ij, W_pot_ij;
kernel_grav_eval(ui, &W_ij);
kernel_grav_force_eval(ui, &W_f_ij);
kernel_grav_pot_eval(ui, &W_pot_ij);
/* Get softened gravity */
*f_ij = mass * h_inv3 * W_ij;
*f_ij = mass * h_inv3 * W_f_ij;
*pot_ij = mass * h_inv * W_pot_ij;
}
/* Get long-range correction */
const float u_lr = r * rlr_inv;
float corr_lr;
kernel_long_grav_eval(u_lr, &corr_lr);
*f_ij *= corr_lr;
float corr_f_lr, corr_pot_lr;
kernel_long_grav_force_eval(u_lr, &corr_f_lr);
kernel_long_grav_pot_eval(u_lr, &corr_pot_lr);
*f_ij *= corr_f_lr;
*pot_ij *= corr_pot_lr;
}
/**
......
src/kernel_gravity.h
View file @ 23129ec0
......@@ -27,14 +27,35 @@
#include "minmax.h"
/**
* @brief Computes the gravity softening function.
* @brief Computes the gravity softening function for potential.
*
* This functions assumes 0 < u < 1.
*
* @param u The ratio of the distance to the softening length $u = x/h$.
* @param W (return) The value of the kernel function $W(x,h)$.
*/
__attribute__((always_inline)) INLINE static void kernel_grav_eval(
__attribute__((always_inline)) INLINE static void kernel_grav_pot_eval(
float u, float *const W) {
/* W(u) = 3u^7 - 15u^6 + 28u^5 - 21u^4 + 7u^2 - 3 */
*W = 3.f * u - 15.f;
*W = *W * u + 28.f;
*W = *W * u - 21.f;
*W = *W * u;
*W = *W * u + 7.f;
*W = *W * u;
*W = *W * u - 3;
}
/**
* @brief Computes the gravity softening function for forces.
*
* This functions assumes 0 < u < 1.
*
* @param u The ratio of the distance to the softening length $u = x/h$.
* @param W (return) The value of the kernel function $W(x,h)$.
*/
__attribute__((always_inline)) INLINE static void kernel_grav_force_eval(
float u, float *const W) {
/* W(u) = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
......
src/kernel_long_gravity.h
View file @ 23129ec0
......@@ -31,12 +31,41 @@
#include <math.h>
/**
* @brief Computes the long-range correction term for the FFT calculation.
* @brief Computes the long-range correction term for the potential calculation
* coming from FFT.
*
* @param u The ratio of the distance to the FFT cell scale \f$u = r/r_s\f$.
* @param W (return) The value of the kernel function.
*/
__attribute__((always_inline)) INLINE static void kernel_long_grav_eval(
__attribute__((always_inline)) INLINE static void kernel_long_grav_pot_eval(
const float u, float *const W) {
#ifdef GADGET2_LONG_RANGE_CORRECTION
const float arg1 = u * 0.5f;
const float term1 = erfcf(arg1);
*W = term1;
#else
const float x = 2.f * u;
const float exp_x = good_approx_expf(x);
const float alpha = 1.f / (1.f + exp_x);
/* We want 2 - 2 exp(x) * alpha */
*W = 1.f - alpha * exp_x;
*W *= 2.f;
#endif
}
/**
* @brief Computes the long-range correction term for the force calculation
* coming from FFT.
*
* @param u The ratio of the distance to the FFT cell scale \f$u = r/r_s\f$.
* @param W (return) The value of the kernel function.
*/
__attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
float u, float *const W) {
#ifdef GADGET2_LONG_RANGE_CORRECTION
......
src/multipole.h
View file @ 23129ec0
......@@ -2384,7 +2384,7 @@ INLINE static void gravity_M2P(const struct multipole *ma,
const float r = r2 * r_inv;
const float u = r * eps_inv;
kernel_grav_eval(u, &W);
kernel_grav_force_eval(u, &W);
/* Get softened gravity */
f = ma->M_000 * eps_inv3 * W;
......
src/runner_doiact_grav.h
View file @ 23129ec0
......@@ -204,8 +204,8 @@ static INLINE void runner_dopair_grav_pp_full(const struct engine *e,
const float h_inv_i = 1.f / h_i;
const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
/* Local accumulators for the acceleration */
float a_x = 0.f, a_y = 0.f, a_z = 0.f;
/* Local accumulators for the acceleration and potential */
float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
/* Make the compiler understand we are in happy vectorization land */
swift_align_information(float, cj_cache->x, SWIFT_CACHE_ALIGNMENT);
......@@ -240,13 +240,15 @@ static INLINE void runner_dopair_grav_pp_full(const struct engine *e,
#endif
/* Interact! */
float f_ij;
runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij);
float f_ij, pot_ij;
runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij,
&pot_ij);
/* Store it back */
a_x -= f_ij * dx;
a_y -= f_ij * dy;
a_z -= f_ij * dz;
pot += pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter if it's not a padded gpart */
......@@ -258,6 +260,7 @@ static INLINE void runner_dopair_grav_pp_full(const struct engine *e,
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
ci_cache->pot[pid] = pot;
}
TIMER_TOC(timer_dopair_grav_pp);
......@@ -295,8 +298,8 @@ static INLINE void runner_dopair_grav_pp_truncated(
const float h_inv_i = 1.f / h_i;
const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
/* Local accumulators for the acceleration */
float a_x = 0.f, a_y = 0.f, a_z = 0.f;
/* Local accumulators for the acceleration and potential */
float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
/* Make the compiler understand we are in happy vectorization land */
swift_align_information(float, cj_cache->x, SWIFT_CACHE_ALIGNMENT);
......@@ -331,14 +334,15 @@ static INLINE void runner_dopair_grav_pp_truncated(
#endif
/* Interact! */
float f_ij;
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
rlr_inv, &f_ij);
rlr_inv, &f_ij, &pot_ij);
/* Store it back */
a_x -= f_ij * dx;
a_y -= f_ij * dy;
a_z -= f_ij * dz;
pot += pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter if it's not a padded gpart */
......@@ -350,6 +354,7 @@ static INLINE void runner_dopair_grav_pp_truncated(
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
ci_cache->pot[pid] = pot;
}
TIMER_TOC(timer_dopair_grav_pp);
......@@ -683,7 +688,7 @@ void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
/* Local accumulators for the acceleration */
float a_x = 0.f, a_y = 0.f, a_z = 0.f;
float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
/* Make the compiler understand we are in happy vectorization land */
swift_align_information(float, ci_cache->x, SWIFT_CACHE_ALIGNMENT);
......@@ -721,13 +726,15 @@ void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
#endif
/* Interact! */
float f_ij;
runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij);
float f_ij, pot_ij;
runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij,
&pot_ij);
/* Store it back */
a_x -= f_ij * dx;
a_y -= f_ij * dy;
a_z -= f_ij * dz;
pot += pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter if it's not a padded gpart */
......@@ -739,6 +746,7 @@ void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
ci_cache->pot[pid] = pot;
}
/* Write back to the particles */
......@@ -808,8 +816,8 @@ void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
const float h_inv_i = 1.f / h_i;
const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
/* Local accumulators for the acceleration */
float a_x = 0.f, a_y = 0.f, a_z = 0.f;
/* Local accumulators for the acceleration and potential */
float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
/* Make the compiler understand we are in happy vectorization land */
swift_align_information(float, ci_cache->x, SWIFT_CACHE_ALIGNMENT);
......@@ -847,14 +855,15 @@ void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
#endif
/* Interact! */
float f_ij;
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
rlr_inv, &f_ij);
rlr_inv, &f_ij, &pot_ij);
/* Store it back */
a_x -= f_ij * dx;
a_y -= f_ij * dy;
a_z -= f_ij * dz;
pot += pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter if it's not a padded gpart */
......@@ -866,6 +875,7 @@ void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
ci_cache->pot[pid] = pot;
}
/* Write back to the particles */
......
Supports Markdown
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