diff --git a/configure.ac b/configure.ac
index 2597d61d01b9df1fa22c1759960179c25b42e226..32ef711f2e306cf16edd6ac5f58a08bd2c172dc1 100644
--- a/configure.ac
+++ b/configure.ac
@@ -917,7 +917,26 @@ case "$with_subgrid" in
;;
esac
+# Gravity scheme.
+AC_ARG_WITH([gravity],
+ [AS_HELP_STRING([--with-gravity=<scheme>],
+ [Gravity scheme to use @<:@default, with-potential, default: default@:>@]
+ )],
+ [with_gravity="$withval"],
+ [with_gravity="default"]
+)
+case "$with_gravity" in
+ with-potential)
+ AC_DEFINE([POTENTIAL_GRAVITY], [1], [Gravity scheme with potential calculation])
+ ;;
+ default)
+ AC_DEFINE([DEFAULT_GRAVITY], [1], [Default gravity scheme])
+ ;;
+ *)
+ AC_MSG_ERROR([Unknown gravity scheme: $with_gravity])
+ ;;
+esac
# Hydro scheme.
AC_ARG_WITH([hydro],
@@ -1351,12 +1370,14 @@ AC_MSG_RESULT([
Equation of state : $with_eos
Adiabatic index : $with_gamma
Riemann solver : $with_riemann
- Cooling function : $with_cooling
- Chemistry : $with_chemistry
- External potential : $with_potential
+ Gravity scheme : $with_gravity
Multipole order : $with_multipole_order
No gravity below ID : $no_gravity_below_id
+ External potential : $with_potential
+
+ Cooling function : $with_cooling
+ Chemistry : $with_chemistry
Individual timers : $enable_timers
Task debugging : $enable_task_debugging
diff --git a/src/debug.c b/src/debug.c
index 05c21de0a73bba3a5e867a4265de0a5c14736a14..da8ef0e118b57a6aa94577898b03bcf7c56b006a 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -62,7 +62,14 @@
#error "Invalid choice of SPH variant"
#endif
+/* Import the right gravity definition */
+#if defined(DEFAULT_GRAVITY)
#include "./gravity/Default/gravity_debug.h"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity_debug.h"
+#else
+#error "Invalid choice of gravity variant"
+#endif
/**
* @brief Looks for the particle with the given id and prints its information to
diff --git a/src/gravity.h b/src/gravity.h
index 6497de8294dfa3f207332ff696ddb992c875eb28..9ff52dc739358e1486e071e268e0f1ab50589e25 100644
--- a/src/gravity.h
+++ b/src/gravity.h
@@ -27,10 +27,19 @@
#include "inline.h"
#include "part.h"
-/* So far only one model here */
-/* Straight-forward import */
+/* Import the right functions */
+#if defined(DEFAULT_GRAVITY)
#include "./gravity/Default/gravity.h"
#include "./gravity/Default/gravity_iact.h"
+#define GRAVITY_IMPLEMENTATION "Default (no potential)"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity.h"
+#include "./gravity/Potential/gravity_iact.h"
+#define GRAVITY_IMPLEMENTATION "With potential calculation"
+#else
+#error "Invalid choice of gravity variant"
+#endif
+
struct engine;
struct space;
diff --git a/src/gravity/Potential/gravity.h b/src/gravity/Potential/gravity.h
new file mode 100644
index 0000000000000000000000000000000000000000..9f0db3f3bde9aef7a847d22dbc36b35b192b9304
--- /dev/null
+++ b/src/gravity/Potential/gravity.h
@@ -0,0 +1,205 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2016 Tom Theuns (tom.theuns@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_DEFAULT_GRAVITY_H
+#define SWIFT_DEFAULT_GRAVITY_H
+
+#include <float.h>
+
+#include "cosmology.h"
+#include "gravity_properties.h"
+#include "kernel_gravity.h"
+#include "minmax.h"
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param gp The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float gravity_get_mass(
+ const struct gpart* restrict gp) {
+
+ return gp->mass;
+}
+
+/**
+ * @brief Returns the softening of a particle
+ *
+ * @param gp The particle of interest
+ * @param grav_props The global gravity properties.
+ */
+__attribute__((always_inline)) INLINE static float gravity_get_softening(
+ const struct gpart* gp, const struct gravity_props* restrict grav_props) {
+
+ return grav_props->epsilon_cur;
+}
+
+/**
+ * @brief Returns the comoving potential of a particle
+ *
+ * @param gp The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_get_comoving_potential(const struct gpart* restrict gp) {
+
+ return gp->potential;
+}
+
+/**
+ * @brief Returns the physical potential of a particle
+ *
+ * @param gp The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_get_physical_potential(const struct gpart* restrict gp,
+ const struct cosmology* cosmo) {
+
+ return gp->potential * cosmo->a_inv;
+}
+
+/**
+ * @brief Computes the gravity time-step of a given particle due to self-gravity
+ *
+ * We use Gadget-2's type 0 time-step criterion.
+ *
+ * @param gp Pointer to the g-particle data.
+ * @param a_hydro The accelerations coming from the hydro scheme (can be 0).
+ * @param grav_props Constants used in the gravity scheme.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_compute_timestep_self(const struct gpart* const gp,
+ const float a_hydro[3],
+ const struct gravity_props* restrict grav_props,
+ const struct cosmology* cosmo) {
+
+ /* Get physical acceleration (gravity contribution) */
+ float a_phys_x = gp->a_grav[0] * cosmo->a_factor_grav_accel;
+ float a_phys_y = gp->a_grav[1] * cosmo->a_factor_grav_accel;
+ float a_phys_z = gp->a_grav[2] * cosmo->a_factor_grav_accel;
+
+ /* Get physical acceleration (hydro contribution) */
+ a_phys_x += a_hydro[0] * cosmo->a_factor_hydro_accel;
+ a_phys_y += a_hydro[1] * cosmo->a_factor_hydro_accel;
+ a_phys_z += a_hydro[2] * cosmo->a_factor_hydro_accel;
+
+ const float ac2 =
+ a_phys_x * a_phys_x + a_phys_y * a_phys_y + a_phys_z * a_phys_z;
+
+ const float ac_inv = (ac2 > 0.f) ? 1.f / sqrtf(ac2) : FLT_MAX;
+
+ const float epsilon = gravity_get_softening(gp, grav_props);
+
+ const float dt = sqrtf(2. * kernel_gravity_softening_plummer_equivalent_inv *
+ cosmo->a * grav_props->eta * epsilon * ac_inv);
+
+ return dt;
+}
+
+/**
+ * @brief Prepares a g-particle for the gravity calculation
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the variaous tasks
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void gravity_init_gpart(
+ struct gpart* gp) {
+
+ /* Zero the acceleration */
+ gp->a_grav[0] = 0.f;
+ gp->a_grav[1] = 0.f;
+ gp->a_grav[2] = 0.f;
+ gp->potential = 0.f;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ gp->potential_PM = 0.f;
+ gp->a_grav_PM[0] = 0.f;
+ gp->a_grav_PM[1] = 0.f;
+ gp->a_grav_PM[2] = 0.f;
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gp->num_interacted = 0;
+#endif
+}
+
+/**
+ * @brief Finishes the gravity calculation.
+ *
+ * Multiplies the forces and accelerations by the appropiate constants
+ *
+ * @param gp The particle to act upon
+ * @param const_G Newton's constant in internal units
+ */
+__attribute__((always_inline)) INLINE static void gravity_end_force(
+ struct gpart* gp, float const_G) {
+
+ /* Let's get physical... */
+ gp->a_grav[0] *= const_G;
+ gp->a_grav[1] *= const_G;
+ gp->a_grav[2] *= const_G;
+ gp->potential *= const_G;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ gp->potential_PM *= const_G;
+ gp->a_grav_PM[0] *= const_G;
+ gp->a_grav_PM[1] *= const_G;
+ gp->a_grav_PM[2] *= const_G;
+#endif
+}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param gp The particle to act upon
+ * @param dt The time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void gravity_kick_extra(
+ struct gpart* gp, float dt) {}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param gp The particle.
+ */
+__attribute__((always_inline)) INLINE static void
+gravity_reset_predicted_values(struct gpart* gp) {}
+
+/**
+ * @brief Initialises the g-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param gp The particle to act upon
+ * @param grav_props The global properties of the gravity calculation.
+ */
+__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
+ struct gpart* gp, const struct gravity_props* grav_props) {
+
+ gp->time_bin = 0;
+
+ gravity_init_gpart(gp);
+}
+
+#endif /* SWIFT_DEFAULT_GRAVITY_H */
diff --git a/src/gravity/Potential/gravity_debug.h b/src/gravity/Potential/gravity_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..dce038c58e1769446861bdf6c9a2a44415642c68
--- /dev/null
+++ b/src/gravity/Potential/gravity_debug.h
@@ -0,0 +1,35 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_DEFAULT_GRAVITY_DEBUG_H
+#define SWIFT_DEFAULT_GRAVITY_DEBUG_H
+
+__attribute__((always_inline)) INLINE static void gravity_debug_particle(
+ const struct gpart* p) {
+ printf(
+ "mass=%.3e time_bin=%d\n"
+ "x=[%.5e,%.5e,%.5e], v_full=[%.5e,%.5e,%.5e], a=[%.5e,%.5e,%.5e]\n",
+ p->mass, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
+ p->v_full[1], p->v_full[2], p->a_grav[0], p->a_grav[1], p->a_grav[2]);
+#ifdef SWIFT_DEBUG_CHECKS
+ printf("num_interacted=%lld ti_drift=%lld ti_kick=%lld\n", p->num_interacted,
+ p->ti_drift, p->ti_kick);
+#endif
+}
+
+#endif /* SWIFT_DEFAULT_GRAVITY_DEBUG_H */
diff --git a/src/gravity/Potential/gravity_iact.h b/src/gravity/Potential/gravity_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..ad477c54f49f9d30ee492d84012d7a9228401c5f
--- /dev/null
+++ b/src/gravity/Potential/gravity_iact.h
@@ -0,0 +1,350 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * 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_DEFAULT_GRAVITY_IACT_H
+#define SWIFT_DEFAULT_GRAVITY_IACT_H
+
+/* Includes. */
+#include "kernel_gravity.h"
+#include "kernel_long_gravity.h"
+#include "multipole.h"
+
+/* Standard headers */
+#include <float.h>
+
+/**
+ * @brief Computes the intensity of the force at a point generated by a
+ * point-mass.
+ *
+ * The returned quantity needs to be multiplied by the distance vector to obtain
+ * the force vector.
+ *
+ * @param r2 Square of the distance to the point-mass.
+ * @param h2 Square of the softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param h_inv3 Cube of the inverse of the softening length.
+ * @param mass Mass of the point-mass.
+ * @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
+ const float r2, const float h2, const float h_inv, const float h_inv3,
+ const float mass, float *f_ij, float *pot_ij) {
+
+ /* Get the inverse distance */
+ const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
+
+ /* Should we soften ? */
+ if (r2 >= h2) {
+
+ /* 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_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_f_ij;
+ *pot_ij = mass * h_inv * W_pot_ij;
+ }
+}
+
+/**
+ * @brief Computes the intensity of the force at a point generated by a
+ * point-mass truncated for long-distance periodicity.
+ *
+ * The returned quantity needs to be multiplied by the distance vector to obtain
+ * the force vector.
+ *
+ * @param r2 Square of the distance to the point-mass.
+ * @param h2 Square of the softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param h_inv3 Cube of the inverse of the softening length.
+ * @param mass Mass of the point-mass.
+ * @param r_s_inv Inverse of the mesh smoothing scale.
+ * @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
+ const float r2, const float h2, const float h_inv, const float h_inv3,
+ const float mass, const float r_s_inv, float *f_ij, float *pot_ij) {
+
+ /* Get the inverse distance */
+ const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
+ const float r = r2 * r_inv;
+
+ /* Should we soften ? */
+ if (r2 >= h2) {
+
+ /* 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_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_f_ij;
+ *pot_ij = mass * h_inv * W_pot_ij;
+ }
+
+ /* Get long-range correction */
+ const float u_lr = r * r_s_inv;
+ 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;
+}
+
+/**
+ * @brief Computes the forces at a point generated by a multipole.
+ *
+ * This assumes M_100 == M_010 == M_001 == 0.
+ * This uses the quadrupole and trace of the octupole terms only and defaults to
+ * the monopole if the code is compiled with low-order gravity only.
+ *
+ * @param r_x x-component of the distance vector to the multipole.
+ * @param r_y y-component of the distance vector to the multipole.
+ * @param r_z z-component of the distance vector to the multipole.
+ * @param r2 Square of the distance vector to the multipole.
+ * @param h The softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param m The multipole.
+ * @param f_x (return) The x-component of the acceleration.
+ * @param f_y (return) The y-component of the acceleration.
+ * @param f_z (return) The z-component of the acceleration.
+ * @param pot (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
+ const float r_x, const float r_y, const float r_z, const float r2,
+ const float h, const float h_inv, const struct multipole *m, float *f_x,
+ float *f_y, float *f_z, float *pot) {
+
+/* In the case where the order is < 3, then there is only a monopole term left.
+ * We can default to the normal P-P interaction with the mass of the multipole
+ * and its CoM as the "particle" property */
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
+
+ float f_ij, pot_ij;
+ runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
+ &f_ij, &pot_ij);
+ *f_x = f_ij * r_x;
+ *f_y = f_ij * r_y;
+ *f_z = f_ij * r_z;
+ *pot = pot_ij;
+
+#else
+
+ /* Get the inverse distance */
+ const float r_inv = 1.f / sqrtf(r2);
+
+ /* Compute the derivatives of the potential */
+ struct potential_derivatives_M2P d;
+ compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, 0.f,
+ &d);
+
+ /* 0th order contributions */
+ *f_x = m->M_000 * d.D_100;
+ *f_y = m->M_000 * d.D_010;
+ *f_z = m->M_000 * d.D_001;
+ *pot = m->M_000 * d.D_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+
+/* 1st order contributions */
+
+/* 1st order contributions are all 0 since the dipole is 0 */
+
+/* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */
+/* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */
+/* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
+/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* 2nd order contributions */
+ *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
+ m->M_101 * d.D_201 + m->M_110 * d.D_210 + m->M_200 * d.D_300;
+ *f_y += m->M_002 * d.D_012 + m->M_011 * d.D_021 + m->M_020 * d.D_030 +
+ m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210;
+ *f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 +
+ m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
+ *pot += m->M_002 * d.D_002 + m->M_011 * d.D_011 + m->M_020 * d.D_020 +
+ m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+
+ /* 3rd order contributions */
+ *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
+ m->M_030 * d.D_130 + m->M_102 * d.D_202 + m->M_111 * d.D_211 +
+ m->M_120 * d.D_220 + m->M_201 * d.D_301 + m->M_210 * d.D_310 +
+ m->M_300 * d.D_400;
+ *f_y += m->M_003 * d.D_013 + m->M_012 * d.D_022 + m->M_021 * d.D_031 +
+ m->M_030 * d.D_040 + m->M_102 * d.D_112 + m->M_111 * d.D_121 +
+ m->M_120 * d.D_130 + m->M_201 * d.D_211 + m->M_210 * d.D_220 +
+ m->M_300 * d.D_310;
+ *f_z += m->M_003 * d.D_004 + m->M_012 * d.D_013 + m->M_021 * d.D_022 +
+ m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 +
+ m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
+ m->M_300 * d.D_301;
+ *pot += m->M_003 * d.D_003 + m->M_012 * d.D_012 + m->M_021 * d.D_021 +
+ m->M_030 * d.D_030 + m->M_102 * d.D_102 + m->M_111 * d.D_111 +
+ m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
+ m->M_300 * d.D_300;
+
+#endif
+
+ /* Take care of the the sign convention */
+ *f_x *= -1.f;
+ *f_y *= -1.f;
+ *f_z *= -1.f;
+ *pot *= -1.f;
+#endif
+}
+
+/**
+ * @brief Computes the forces at a point generated by a multipole, truncated for
+ * long-range periodicity.
+ *
+ * This assumes M_100 == M_010 == M_001 == 0.
+ * This uses the quadrupole term and trace of the octupole terms only and
+ * defaults to the monopole if the code is compiled with low-order gravity only.
+ *
+ * @param r_x x-component of the distance vector to the multipole.
+ * @param r_y y-component of the distance vector to the multipole.
+ * @param r_z z-component of the distance vector to the multipole.
+ * @param r2 Square of the distance vector to the multipole.
+ * @param h The softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param r_s_inv The inverse of the gravity mesh-smoothing scale.
+ * @param m The multipole.
+ * @param f_x (return) The x-component of the acceleration.
+ * @param f_y (return) The y-component of the acceleration.
+ * @param f_z (return) The z-component of the acceleration.
+ * @param pot (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
+ const float r_x, const float r_y, const float r_z, const float r2,
+ const float h, const float h_inv, const float r_s_inv,
+ const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) {
+
+/* In the case where the order is < 3, then there is only a monopole term left.
+ * We can default to the normal P-P interaction with the mass of the multipole
+ * and its CoM as the "particle" property */
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
+
+ float f_ij, pot_ij;
+ runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
+ m->M_000, r_s_inv, &f_ij, &pot_ij);
+ *f_x = f_ij * r_x;
+ *f_y = f_ij * r_y;
+ *f_z = f_ij * r_z;
+ *pot = -pot_ij;
+
+#else
+
+ /* Get the inverse distance */
+ const float r_inv = 1.f / sqrtf(r2);
+
+ /* Compute the derivatives of the potential */
+ struct potential_derivatives_M2P d;
+ compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
+ r_s_inv, &d);
+
+ /* 0th order contributions */
+ *f_x = m->M_000 * d.D_100;
+ *f_y = m->M_000 * d.D_010;
+ *f_z = m->M_000 * d.D_001;
+ *pot = m->M_000 * d.D_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+
+/* 1st order contributions */
+
+/* 1st order contributions are all 0 since the dipole is 0 */
+
+/* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */
+/* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */
+/* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
+/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* 2nd order contributions */
+ *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
+ m->M_101 * d.D_201 + m->M_110 * d.D_210 + m->M_200 * d.D_300;
+ *f_y += m->M_002 * d.D_012 + m->M_011 * d.D_021 + m->M_020 * d.D_030 +
+ m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210;
+ *f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 +
+ m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
+ *pot += m->M_002 * d.D_002 + m->M_011 * d.D_011 + m->M_020 * d.D_020 +
+ m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
+
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+
+ /* 3rd order contributions */
+ *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
+ m->M_030 * d.D_130 + m->M_102 * d.D_202 + m->M_111 * d.D_211 +
+ m->M_120 * d.D_220 + m->M_201 * d.D_301 + m->M_210 * d.D_310 +
+ m->M_300 * d.D_400;
+ *f_y += m->M_003 * d.D_013 + m->M_012 * d.D_022 + m->M_021 * d.D_031 +
+ m->M_030 * d.D_040 + m->M_102 * d.D_112 + m->M_111 * d.D_121 +
+ m->M_120 * d.D_130 + m->M_201 * d.D_211 + m->M_210 * d.D_220 +
+ m->M_300 * d.D_310;
+ *f_z += m->M_003 * d.D_004 + m->M_012 * d.D_013 + m->M_021 * d.D_022 +
+ m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 +
+ m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
+ m->M_300 * d.D_301;
+ *pot += m->M_003 * d.D_003 + m->M_012 * d.D_012 + m->M_021 * d.D_021 +
+ m->M_030 * d.D_030 + m->M_102 * d.D_102 + m->M_111 * d.D_111 +
+ m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
+ m->M_300 * d.D_300;
+
+#endif
+
+ /* Take care of the the sign convention */
+ *f_x *= -1.f;
+ *f_y *= -1.f;
+ *f_z *= -1.f;
+ *pot *= -1.f;
+#endif
+}
+
+#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
diff --git a/src/gravity/Potential/gravity_io.h b/src/gravity/Potential/gravity_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..7f453179641e2ba16b30e3172ddd7853245a1d2f
--- /dev/null
+++ b/src/gravity/Potential/gravity_io.h
@@ -0,0 +1,122 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_DEFAULT_GRAVITY_IO_H
+#define SWIFT_DEFAULT_GRAVITY_IO_H
+
+#include "io_properties.h"
+
+INLINE static void convert_gpart_pos(const struct engine* e,
+ const struct gpart* gp, double* ret) {
+
+ if (e->s->periodic) {
+ ret[0] = box_wrap(gp->x[0], 0.0, e->s->dim[0]);
+ ret[1] = box_wrap(gp->x[1], 0.0, e->s->dim[1]);
+ ret[2] = box_wrap(gp->x[2], 0.0, e->s->dim[2]);
+ } else {
+ ret[0] = gp->x[0];
+ ret[1] = gp->x[1];
+ ret[2] = gp->x[2];
+ }
+}
+
+INLINE static void convert_gpart_vel(const struct engine* e,
+ const struct gpart* gp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, gp->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ ret[0] = gp->v_full[0] + gp->a_grav[0] * dt_kick_grav;
+ ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav;
+ ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav;
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a_inv;
+ ret[1] *= cosmo->a_inv;
+ ret[2] *= cosmo->a_inv;
+}
+
+/**
+ * @brief Specifies which g-particle fields to read from a dataset
+ *
+ * @param gparts The g-particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+INLINE static void darkmatter_read_particles(struct gpart* gparts,
+ struct io_props* list,
+ int* num_fields) {
+
+ /* Say how much we want to read */
+ *num_fields = 4;
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ UNIT_CONV_LENGTH, gparts, x);
+ list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+ UNIT_CONV_SPEED, gparts, v_full);
+ list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ gparts, mass);
+ list[3] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset);
+}
+
+/**
+ * @brief Specifies which g-particle fields to write to a dataset
+ *
+ * @param gparts The g-particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+INLINE static void darkmatter_write_particles(const struct gpart* gparts,
+ struct io_props* list,
+ int* num_fields) {
+
+ /* Say how much we want to write */
+ *num_fields = 5;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field_convert_gpart(
+ "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, gparts, convert_gpart_pos);
+ list[1] = io_make_output_field_convert_gpart(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, gparts, convert_gpart_vel);
+ list[2] =
+ io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, gparts, mass);
+ list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+ UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset);
+ list[4] = io_make_output_field("Potential", FLOAT, 1, UNIT_CONV_POTENTIAL,
+ gparts, potential);
+}
+
+#endif /* SWIFT_DEFAULT_GRAVITY_IO_H */
diff --git a/src/gravity/Potential/gravity_part.h b/src/gravity/Potential/gravity_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..d325baf4de0938d8df539d38ae10f8f3f3ec7d2b
--- /dev/null
+++ b/src/gravity/Potential/gravity_part.h
@@ -0,0 +1,83 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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_DEFAULT_GRAVITY_PART_H
+#define SWIFT_DEFAULT_GRAVITY_PART_H
+
+/* Gravity particle. */
+struct gpart {
+
+ /*! Particle ID. If negative, it is the negative offset of the #part with
+ which this gpart is linked. */
+ long long id_or_neg_offset;
+
+ /*! Particle position. */
+ double x[3];
+
+ /*! Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
+
+ /*! Particle velocity. */
+ float v_full[3];
+
+ /*! Particle acceleration. */
+ float a_grav[3];
+
+ /*! Particle mass. */
+ float mass;
+
+ /*! Gravitational potential */
+ float potential;
+
+ /*! Time-step length */
+ timebin_t time_bin;
+
+ /*! Type of the #gpart (DM, gas, star, ...) */
+ enum part_type type;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Numer of gparts this gpart interacted with */
+ long long num_interacted;
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+#endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+ /*! Acceleration taken from the mesh only */
+ float a_grav_PM[3];
+
+ /*! Potential taken from the mesh only */
+ float potential_PM;
+
+ /* Brute-force particle acceleration. */
+ double a_grav_exact[3];
+
+ /* Brute-force particle potential. */
+ double potential_exact;
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_DEFAULT_GRAVITY_PART_H */
diff --git a/src/gravity_io.h b/src/gravity_io.h
index 6276e50473de64abd1014bdb36a63a14e02ca8cf..752ee906081d706865e62bae7c8f505e9ca64347 100644
--- a/src/gravity_io.h
+++ b/src/gravity_io.h
@@ -19,8 +19,19 @@
#ifndef SWIFT_GRAVITY_IO_H
#define SWIFT_GRAVITY_IO_H
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
#include "./const.h"
+/* Import the right functions */
+#if defined(DEFAULT_GRAVITY)
#include "./gravity/Default/gravity_io.h"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity_io.h"
+#else
+#error "Invalid choice of gravity variant"
+#endif
#endif /* SWIFT_GRAVITY_IO_H */
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 22856c25fbbaebd1dd74c5592ce8fd914e76b61c..fc1ce1d62e02c32d44667d602448fc4eb3a65344 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -112,6 +112,8 @@ void gravity_update(struct gravity_props *p, const struct cosmology *cosmo) {
void gravity_props_print(const struct gravity_props *p) {
+ message("Self-gravity scheme: %s", GRAVITY_IMPLEMENTATION);
+
message("Self-gravity scheme: FMM-MM with m-poles of order %d",
SELF_GRAVITY_MULTIPOLE_ORDER);
@@ -167,6 +169,7 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
"Maximal physical softening length (Plummer equivalent)",
p->epsilon_max_physical);
io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
+ io_write_attribute_s(h_grpgrav, "Scheme", GRAVITY_IMPLEMENTATION);
io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
io_write_attribute_f(h_grpgrav, "Mesh r_cut_max ratio", p->r_cut_max_ratio);
diff --git a/src/part.h b/src/part.h
index 03ec331cb17b95b0133be568d6e857a44d1eaf73..bca84cc0212e79e15ffbeeeb0bbcfc714d5481be 100644
--- a/src/part.h
+++ b/src/part.h
@@ -77,7 +77,13 @@
#endif
/* Import the right gravity particle definition */
+#if defined(DEFAULT_GRAVITY)
#include "./gravity/Default/gravity_part.h"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity_part.h"
+#else
+#error "Invalid choice of gravity variant"
+#endif
/* Import the right star particle definition */
#include "./stars/Default/star_part.h"