diff --git a/src/const.h b/src/const.h
index 8914febf7a7aecddbff65386ae2ceef9da060db1..09ec7636d1bcedc50af21d13ebd34f6eff940ee8 100644
--- a/src/const.h
+++ b/src/const.h
@@ -41,9 +41,6 @@
/* Self gravity stuff. */
#define const_gravity_multipole_order 1
-#define const_gravity_a_smooth FLT_MAX
-#define const_gravity_r_cut 4.5f
-#define const_gravity_eta 0.025f
/* Type of gradients to use (GIZMO_SPH only) */
/* If no option is chosen, no gradients are used (first order scheme) */
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 048bc3ee832d4cf9f939c41a0056ca490729f9fa..93a65e2f5a70e09ad14280cf9c334753359fb8b5 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -21,15 +21,18 @@
#define SWIFT_DEFAULT_GRAVITY_H
#include <float.h>
+#include "gravity_properties.h"
#include "minmax.h"
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
* @param gp Pointer to the g-particle data.
+ * @param grav_props Constants used in the gravity scheme.
*/
__attribute__((always_inline)) INLINE static float
-gravity_compute_timestep_self(const struct gpart* const gp) {
+gravity_compute_timestep_self(const struct gpart* const gp,
+ const struct gravity_props* restrict grav_props) {
const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
gp->a_grav[1] * gp->a_grav[1] +
@@ -37,7 +40,7 @@ gravity_compute_timestep_self(const struct gpart* const gp) {
const float ac = (ac2 > 0.f) ? sqrtf(ac2) : FLT_MIN;
- const float dt = sqrtf(2.f * const_gravity_eta * gp->epsilon / ac);
+ const float dt = sqrtf(2.f * grav_props->eta * gp->epsilon / ac);
return dt;
}
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 9498caeb0b4b15cc849160375a50ef984393d847..4f2aa3f5868153644b21b7d1b971103b00ac6d25 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -69,7 +69,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
const int gcount = ci->gcount;
struct gpart *restrict gparts = ci->gparts;
const struct multipole *multi = cj->multipole;
- const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
+ const float a_smooth = e->gravity_properties->a_smooth;
+ const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
TIMER_TIC;
@@ -138,7 +139,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
const int gcount_j = cj->gcount;
struct gpart *restrict gparts_i = ci->gparts;
struct gpart *restrict gparts_j = cj->gparts;
- const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
+ const float a_smooth = e->gravity_properties->a_smooth;
+ const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
TIMER_TIC;
@@ -244,7 +246,8 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
const struct engine *e = r->e;
const int gcount = c->gcount;
struct gpart *restrict gparts = c->gparts;
- const float rlr_inv = 1. / (const_gravity_a_smooth * c->super->width[0]);
+ const float a_smooth = e->gravity_properties->a_smooth;
+ const float rlr_inv = 1. / (a_smooth * c->super->width[0]);
TIMER_TIC;
diff --git a/src/timestep.h b/src/timestep.h
index 6497e25984ef36d759afd12616b45c8166816dfb..99ca1b4f90cc7b8894d4dfe0e0d2670da8f01d65 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -78,7 +78,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_gpart_timestep(
e->physical_constants, gp));
if (e->policy & engine_policy_self_gravity)
- new_dt = min(new_dt, gravity_compute_timestep_self(gp));
+ new_dt =
+ min(new_dt, gravity_compute_timestep_self(gp, e->gravity_properties));
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
@@ -121,7 +122,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
e->physical_constants, p->gpart));
if (e->policy & engine_policy_self_gravity)
- new_dt_grav = min(new_dt_grav, gravity_compute_timestep_self(p->gpart));
+ new_dt_grav = min(new_dt_grav, gravity_compute_timestep_self(
+ p->gpart, e->gravity_properties));
}
/* Final time-step is minimum of hydro and gravity */
@@ -163,7 +165,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
e->physical_constants, sp->gpart));
if (e->policy & engine_policy_self_gravity)
- new_dt = min(new_dt, gravity_compute_timestep_self(sp->gpart));
+ new_dt = min(new_dt, gravity_compute_timestep_self(sp->gpart,
+ e->gravity_properties));
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
diff --git a/src/tools.c b/src/tools.c
index ab11d1f5930cf5319aaf6424f1559f144718e154..89ac286fb435c01b361bdea66e62dd2d7f41ee24 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -732,19 +732,24 @@ int compare_particles(struct part a, struct part b, double threshold) {
#endif
}
-/** @brief Computes the forces between all g-particles using the N^2 algorithm
+/**
+ * @brief Computes the forces between all g-particles using the N^2 algorithm
*
* Overwrites the accelerations of the gparts with the values.
* Do not use for actual runs.
*
* @brief gparts The array of particles.
* @brief gcount The number of particles.
+ * @brief constants Physical constants in internal units.
+ * @brief gravity_properties Constants governing the gravity scheme.
*/
void gravity_n2(struct gpart *gparts, const int gcount,
- const struct phys_const *constants, float rlr) {
+ const struct phys_const *constants,
+ const struct gravity_props *gravity_properties, float rlr) {
const float rlr_inv = 1. / rlr;
- const float max_d = const_gravity_r_cut * rlr;
+ const float r_cut = gravity_properties->r_cut;
+ const float max_d = r_cut * rlr;
const float max_d2 = max_d * max_d;
message("rlr_inv= %f", rlr_inv);
diff --git a/src/tools.h b/src/tools.h
index ece3078dce7cc8ab4b15538a1e5d9a990d81b36d..4d9e8d3ef86f9ad2661118acf008797893ea5bd7 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -23,6 +23,7 @@
#define SWIFT_TOOL_H
#include "cell.h"
+#include "gravity_properties.h"
#include "part.h"
#include "physical_constants.h"
#include "runner.h"
@@ -45,8 +46,8 @@ void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic);
double random_uniform(double a, double b);
void shuffle_particles(struct part *parts, const int count);
void gravity_n2(struct gpart *gparts, const int gcount,
- const struct phys_const *constants, float rlr);
-
+ const struct phys_const *constants,
+ const struct gravity_props *gravity_properties, float rlr);
int compare_values(double a, double b, double threshold, double *absDiff,
double *absSum, double *relDiff);
int compare_particles(struct part a, struct part b, double threshold);