diff --git a/src/Makefile.am b/src/Makefile.am
index 9b043c829a20e99306eac7fa7b352859248ab7ae..fe17d77a874f3b2780b1ae00ba7523c3e63f51e6 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -64,12 +64,14 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
- kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h \
+ gravity_iact.h kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h \
runner_doiact_nosort.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
dimension.h part_type.h periodic.h memswap.h dump.h logger.h sign.h \
gravity.h gravity_io.h gravity_cache.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 \
+ gravity/Potential/gravity.h gravity/Potential/gravity_iact.h gravity/Potential/gravity_io.h \
+ gravity/Potential/gravity_debug.h gravity/Potential/gravity_part.h \
sourceterms.h \
equation_of_state.h \
equation_of_state/ideal_gas/equation_of_state.h equation_of_state/isothermal/equation_of_state.h \
diff --git a/src/gravity.h b/src/gravity.h
index 9ff52dc739358e1486e071e268e0f1ab50589e25..57cad8eba5f0772f85affd031a450c3ecb4dde59 100644
--- a/src/gravity.h
+++ b/src/gravity.h
@@ -30,17 +30,14 @@
/* 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/Default/gravity.h b/src/gravity/Default/gravity.h
index 9f0db3f3bde9aef7a847d22dbc36b35b192b9304..c90ca39e53c525635fe328312babf84f266dfd9a 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -51,19 +51,36 @@ __attribute__((always_inline)) INLINE static float gravity_get_softening(
}
/**
- * @brief Returns the comoving potential of a particle
+ * @brief Add a contribution to this particle's potential.
+ *
+ * Here we do nothing as this version does not accumulate potential.
+ *
+ * @param gp The particle.
+ * @param pot The contribution to add.
+ */
+__attribute__((always_inline)) INLINE static void
+gravity_add_comoving_potential(struct gpart* restrict gp, float pot) {}
+
+/**
+ * @brief Returns the comoving potential of a particle.
+ *
+ * This returns 0 as this flavour of gravity does not compute the
+ * particles' potential.
*
* @param gp The particle of interest
*/
__attribute__((always_inline)) INLINE static float
gravity_get_comoving_potential(const struct gpart* restrict gp) {
- return gp->potential;
+ return 0.f;
}
/**
* @brief Returns the physical potential of a particle
*
+ * This returns 0 as this flavour of gravity does not compute the
+ * particles' potential.
+ *
* @param gp The particle of interest.
* @param cosmo The cosmological model.
*/
@@ -71,7 +88,7 @@ __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;
+ return 0.f;
}
/**
@@ -128,7 +145,6 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
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;
@@ -157,7 +173,6 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
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;
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index ad477c54f49f9d30ee492d84012d7a9228401c5f..6fd7ae205ecb124344f0780a5e8a276445ad07b4 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -55,21 +55,21 @@ __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_f_ij, W_pot_ij;
+ float W_f_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;
}
+
+ /* No potential calculation */
+ *pot_ij = 0.f;
}
/**
@@ -101,28 +101,26 @@ __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_f_ij, W_pot_ij;
+ float W_f_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;
+ float corr_f_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;
+
+ /* No potential calculation */
+ *pot_ij = 0.f;
}
/**
@@ -160,7 +158,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
*f_x = f_ij * r_x;
*f_y = f_ij * r_y;
*f_z = f_ij * r_z;
- *pot = pot_ij;
#else
@@ -176,7 +173,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
*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
@@ -187,7 +183,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
/* *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
@@ -200,8 +195,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
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
@@ -220,10 +213,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
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
@@ -231,8 +220,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
*f_x *= -1.f;
*f_y *= -1.f;
*f_z *= -1.f;
- *pot *= -1.f;
#endif
+
+ /* No potential calculation */
+ *pot = 0.f;
}
/**
@@ -272,7 +263,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
*f_x = f_ij * r_x;
*f_y = f_ij * r_y;
*f_z = f_ij * r_z;
- *pot = -pot_ij;
#else
@@ -288,7 +278,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
*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
@@ -299,7 +288,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
/* *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
@@ -312,8 +300,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
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
@@ -332,19 +318,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
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
+
+ /* No potential calculation */
+ *pot = 0.f;
}
#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index 7f453179641e2ba16b30e3172ddd7853245a1d2f..2150b3c445eb8fbd109e209508269e81d2e773c3 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/src/gravity/Default/gravity_io.h
@@ -115,8 +115,6 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts,
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/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index d325baf4de0938d8df539d38ae10f8f3f3ec7d2b..9af294715bb8ae9db4486a347fbde3e72f55fe82 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -41,9 +41,6 @@ struct gpart {
/*! Particle mass. */
float mass;
- /*! Gravitational potential */
- float potential;
-
/*! Time-step length */
timebin_t time_bin;
diff --git a/src/gravity/Potential/gravity.h b/src/gravity/Potential/gravity.h
index 9f0db3f3bde9aef7a847d22dbc36b35b192b9304..e536a221679be2553ebda54d46419cc275371620 100644
--- a/src/gravity/Potential/gravity.h
+++ b/src/gravity/Potential/gravity.h
@@ -50,6 +50,18 @@ __attribute__((always_inline)) INLINE static float gravity_get_softening(
return grav_props->epsilon_cur;
}
+/**
+ * @brief Add a contribution to this particle's potential.
+ *
+ * @param gp The particle.
+ * @param pot The contribution to add.
+ */
+__attribute__((always_inline)) INLINE static void
+gravity_add_comoving_potential(struct gpart* restrict gp, float pot) {
+
+ gp->potential += pot;
+}
+
/**
* @brief Returns the comoving potential of a particle
*
diff --git a/src/gravity_cache.h b/src/gravity_cache.h
index 9101fc763d69ace10c57536ebf59f655744b212b..821f044429b445c28ff8ae39b8dc65304dd2b42d 100644
--- a/src/gravity_cache.h
+++ b/src/gravity_cache.h
@@ -429,7 +429,7 @@ gravity_cache_populate_all_mpole(const timebin_t max_active_bin,
* @param gparts The #gpart array to write to.
* @param gcount The number of particles to write.
*/
-__attribute__((always_inline)) INLINE void gravity_cache_write_back(
+__attribute__((always_inline)) INLINE static void gravity_cache_write_back(
const struct gravity_cache *c, struct gpart *restrict gparts,
const int gcount) {
@@ -446,7 +446,7 @@ __attribute__((always_inline)) INLINE void gravity_cache_write_back(
gparts[i].a_grav[0] += a_x[i];
gparts[i].a_grav[1] += a_y[i];
gparts[i].a_grav[2] += a_z[i];
- gparts[i].potential += pot[i];
+ gravity_add_comoving_potential(&gparts[i], pot[i]);
}
}
}
diff --git a/src/gravity_iact.h b/src/gravity_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..2bd2272128ad430a65ab723847481f5684a30a71
--- /dev/null
+++ b/src/gravity_iact.h
@@ -0,0 +1,40 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 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_GRAVITY_IACT_H
+#define SWIFT_GRAVITY_IACT_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+#include "inline.h"
+#include "part.h"
+
+/* Import the right functions */
+#if defined(DEFAULT_GRAVITY)
+#include "./gravity/Default/gravity_iact.h"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity_iact.h"
+#else
+#error "Invalid choice of gravity variant"
+#endif
+
+
+#endif
diff --git a/src/mesh_gravity.c b/src/mesh_gravity.c
index 49e66763f7f9ef47ca79d0024fa65b9ce7803fc1..1ac0f8ea6b3b414c1338573cb5a50c3d0f3eb1b9 100644
--- a/src/mesh_gravity.c
+++ b/src/mesh_gravity.c
@@ -254,7 +254,7 @@ void mesh_to_gparts_CIC(struct gpart* gp, const double* pot, int N, double fac,
/* ---- */
/* Store things back */
- gp->potential += p;
+ gravity_add_comoving_potential(gp, p);
gp->a_grav[0] += fac * a[0];
gp->a_grav[1] += fac * a[1];
gp->a_grav[2] += fac * a[2];
diff --git a/src/multipole.h b/src/multipole.h
index ac746728ee39d597c92bba8c260140c10adbddd1..51fa99dde3c54f8c3190bf564e934e2de1b3d1e4 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -31,6 +31,7 @@
#include "align.h"
#include "const.h"
#include "error.h"
+#include "gravity.h"
#include "gravity_derivatives.h"
#include "gravity_properties.h"
#include "gravity_softened_derivatives.h"
@@ -2380,7 +2381,7 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
gp->a_grav[0] += a_grav[0];
gp->a_grav[1] += a_grav[1];
gp->a_grav[2] += a_grav[2];
- gp->potential += pot;
+ gravity_add_comoving_potential(gp, pot);
}
/**
diff --git a/src/runner.c b/src/runner.c
index fe8e7483e24406bf23bfec6a684814c664158a69..f96c113b0d906263bd9a2ef10a1354df24a0a585 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1654,7 +1654,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
- const struct space *s = e->s;
+ //const struct space *s = e->s;
const struct cosmology *cosmo = e->cosmology;
const int count = c->count;
const int gcount = c->gcount;
@@ -1662,12 +1662,12 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct gpart *restrict gparts = c->gparts;
struct spart *restrict sparts = c->sparts;
- const int periodic = s->periodic;
- const int with_cosmology = e->policy & engine_policy_cosmology;
- const float Omega_m = e->cosmology->Omega_m;
- const float H0 = e->cosmology->H0;
+ //const int periodic = s->periodic;
+ //const int with_cosmology = e->policy & engine_policy_cosmology;
+ //const float Omega_m = e->cosmology->Omega_m;
+ //const float H0 = e->cosmology->H0;
const float G_newton = e->physical_constants->const_newton_G;
- const float rho_crit0 = 3.f * H0 * H0 / (8.f * M_PI * G_newton);
+ //const float rho_crit0 = 3.f * H0 * H0 / (8.f * M_PI * G_newton);
TIMER_TIC;
@@ -1704,16 +1704,17 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
/* Finish the force calculation */
gravity_end_force(gp, G_newton);
- /* Apply periodic BC contribution to the potential */
- if (with_cosmology && periodic) {
- const float mass = gravity_get_mass(gp);
- const float mass2 = mass * mass;
-
- /* This correction term matches the one used in Gadget-2 */
- /* The numerical constant is taken from Hernquist, Bouchet & Suto 1991
- */
- gp->potential -= 2.8372975f * cbrtf(mass2 * Omega_m * rho_crit0);
- }
+ // MATTHIEU: Move this and change it!
+ /* /\* Apply periodic BC contribution to the potential *\/ */
+ /* if (with_cosmology && periodic) { */
+ /* const float mass = gravity_get_mass(gp); */
+ /* const float mass2 = mass * mass; */
+
+ /* /\* This correction term matches the one used in Gadget-2 *\/ */
+ /* /\* The numerical constant is taken from Hernquist, Bouchet & Suto 1991 */
+ /* *\/ */
+ /* gp->potential -= 2.8372975f * cbrtf(mass2 * Omega_m * rho_crit0); */
+ /* } */
#ifdef SWIFT_NO_GRAVITY_BELOW_ID
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index bc0599df9ae8e311ab571d95529232d951860323..654640ef8b99e51f206af9f50ed9aa89cd1e0bd2 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -24,6 +24,7 @@
#include "active.h"
#include "cell.h"
#include "gravity.h"
+#include "gravity_iact.h"
#include "gravity_cache.h"
#include "inline.h"
#include "part.h"
@@ -1010,15 +1011,6 @@ static INLINE void runner_doself_grav_pp_truncated(
runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
r_s_inv, &f_ij, &pot_ij);
- /* if (e->s->parts[-gparts[pid].id_or_neg_offset].id == ICHECK) { */
- /* if (pjd < gcount) */
- /* message("Interacting with particle ID= %lld f_ij=%e", */
- /* e->s->parts[-gparts[pjd].id_or_neg_offset].id, f_ij); */
- /* // else */
- /* // message("Interacting with particle ID= %lld (padded) f_ij=%e", */
- /* // e->s->parts[-gparts[pjd].id_or_neg_offset].id, f_ij); */
- /* } */
-
/* Store it back */
a_x += f_ij * dx;
a_y += f_ij * dy;