Do not compute the potential of the particles in the Default gravity model.

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 \

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;
 

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;

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 */

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 */

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;
 

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
  *

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]);
     }
   }
 }

src/gravity_iact.h

+/*******************************************************************************
+ * 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

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];

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);
 }
 
 /**

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
 

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;