diff --git a/src/Makefile.am b/src/Makefile.am index cd5a47ea8df4dc4826de8849b3692d9d8a4c53f6..311c49570a06ff8c269cc26bc8986de291e86166 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -47,8 +47,8 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ # Include files for distribution, not installation. nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \ - vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h kick.h \ - timestep.h drift.h \ + kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h \ + minmax.h kick.h timestep.h drift.h \ gravity.h gravity_io.h \ gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \ gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \ diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h index d0624ab2b23f6098aa21f32dda81bc83e6f2489c..997bff0fb1835dd2574ae7ef91e8b283d38e8450 100644 --- a/src/gravity/Default/gravity_iact.h +++ b/src/gravity/Default/gravity_iact.h @@ -23,6 +23,7 @@ /* Includes. */ #include "const.h" #include "kernel_gravity.h" +#include "kernel_long_gravity.h" #include "multipole.h" #include "vector.h" @@ -49,7 +50,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp( if (r >= hi) { - /* Get Newtonian graavity */ + /* Get Newtonian gravity */ fi = mj * ir * ir * ir; } else { @@ -61,7 +62,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp( if (r >= hj) { - /* Get Newtonian graavity */ + /* Get Newtonian gravity */ fj = mi * ir * ir * ir; } else { @@ -102,7 +103,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym( if (r >= hi) { - /* Get Newtonian graavity */ + /* Get Newtonian gravity */ f = mj * ir * ir * ir; } else { diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h index 1e39f96854fd5b3c325edaa2e597a1a6f3b81f60..b38feb5758debf87add2007ee3684d869f393f7e 100644 --- a/src/kernel_gravity.h +++ b/src/kernel_gravity.h @@ -63,9 +63,9 @@ static const float 0.f}; /* 1 < u */ /** - * @brief Computes the kernel function. + * @brief Computes the gravity softening function. * - * @param u The ratio of the distance to the smoothing length $u = x/h$. + * @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( diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h new file mode 100644 index 0000000000000000000000000000000000000000..d247c7a461d4bd116f30ab106143f6c75e1b941e --- /dev/null +++ b/src/kernel_long_gravity.h @@ -0,0 +1,50 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * + * 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_KERNEL_LONG_GRAVITY_H +#define SWIFT_KERNEL_LONG_GRAVITY_H + +#include <math.h> + +/* Includes. */ +#include "const.h" +#include "inline.h" +#include "vector.h" + +#define one_over_sqrt_pi ((float)(M_2_SQRTPI * 0.5)) + +/** + * @brief Computes the long-range correction term for the FFT calculation. + * + * @param u The ratio of the distance to the FFT cell scale $u = x/A$. + * @param W (return) The value of the kernel function. + */ +__attribute__((always_inline)) INLINE static void kernel_long_grav_eval( + float u, float *const W) { + + const float arg1 = u * 0.5f; + const float arg2 = u * one_over_sqrt_pi; + const float arg3 = -arg1 * arg1; + + const float term1 = erfc(arg1); + const float term2 = arg2 * expf(arg3); + + *W = term1 + term2; +} + +#endif // SWIFT_KERNEL_LONG_GRAVITY_H