diff --git a/src/const.h b/src/const.h
index 6b9a4d2f723a127caa1eff90d1b27b2d1acf6915..c004e4e2e6438d578c63fef04486dbee2c418eb1 100644
--- a/src/const.h
+++ b/src/const.h
@@ -41,9 +41,7 @@
/* Gravity stuff. */
#define multipole_order 2
-#define const_theta_max 0.57735f
-#define const_G 6.672e-8f /* Gravitational constant. */
-#define const_epsilon 0.0014f /* Gravity blending distance. */
+#define const_gravity_eta 0.025f
/* Kernel function to use */
#define CUBIC_SPLINE_KERNEL
diff --git a/src/engine.c b/src/engine.c
index 5b7bef4121f23eb871899a0e784b40248735da12..de78f1ba34776a2b9987674801ed35f3c5498f1f 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2262,7 +2262,7 @@ void engine_step(struct engine *e) {
FILE *file = fopen("grav_swift.dat", "w");
for (size_t k = 0; k < s->nr_gparts; ++k) {
- fprintf(file, "%lld %f %f %f %f %f %f\n", s->gparts[k].id,
+ fprintf(file, "%lld %f %f %f %e %e %e\n", s->gparts[k].id,
s->gparts[k].x[0], s->gparts[k].x[1], s->gparts[k].x[2],
s->gparts[k].a_grav[0], s->gparts[k].a_grav[1],
s->gparts[k].a_grav[2]);
@@ -2275,10 +2275,10 @@ void engine_step(struct engine *e) {
/* Check the gravity accelerations */
struct gpart *temp = malloc(s->nr_gparts * sizeof(struct gpart));
memcpy(temp, s->gparts, s->nr_gparts * sizeof(struct gpart));
- gravity_n2(temp, s->nr_gparts);
+ gravity_n2(temp, s->nr_gparts, e->physical_constants);
file = fopen("grav_brute.dat", "w");
for (size_t k = 0; k < s->nr_gparts; ++k) {
- fprintf(file, "%lld %f %f %f %f %f %f\n", temp[k].id, temp[k].x[0],
+ fprintf(file, "%lld %f %f %f %e %e %e\n", temp[k].id, temp[k].x[0],
temp[k].x[1], temp[k].x[2], temp[k].a_grav[0], temp[k].a_grav[1],
temp[k].a_grav[2]);
}
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index cd3e0dc80048cd8f74e37eb347e4cc55b5d2c0c3..e5e40def18a50b89e9ec1eece150536fdf939fc3 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -50,8 +50,6 @@ __attribute__((always_inline))
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
- * This function only branches towards the potential chosen by the user.
- *
* @param phys_const The physical constants in internal units.
* @param gp Pointer to the g-particle data.
*/
@@ -60,7 +58,13 @@ __attribute__((always_inline))
const struct phys_const* const phys_const,
const struct gpart* const gp) {
- float dt = FLT_MAX;
+ const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
+ gp->a_grav[1] * gp->a_grav[1] +
+ gp->a_grav[2] * gp->a_grav[2];
+
+ const float ac = ac2 > 0. ? sqrtf(ac) : FLT_MIN;
+
+ const float dt = sqrt(2.f * const_gravity_eta * gp->epsilon / ac);
return dt;
}
@@ -101,9 +105,16 @@ __attribute__((always_inline))
* Multiplies the forces and accelerations by the appropiate constants
*
* @param gp The particle to act upon
+ * @param const_G Newton's constant
*/
__attribute__((always_inline))
- INLINE static void gravity_end_force(struct gpart* gp) {}
+INLINE static void gravity_end_force(struct gpart* gp, double const_G) {
+
+ /* Let's get physical... */
+ gp->a_grav[0] *= const_G;
+ gp->a_grav[1] *= const_G;
+ gp->a_grav[2] *= const_G;
+}
/**
* @brief Computes the gravitational acceleration induced by external potentials
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 62023345f174eb8cb9bae4d4438bdd50c9969494..6d8aa3261fe30b765a58755856f29c98b0c10fd8 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -29,88 +29,12 @@
* @brief Gravity potential
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav(
- float r2, float *dx, struct gpart *pi, struct gpart *pj) {
-
- float ir, r;
- float w, acc;
- float mi = pi->mass, mj = pj->mass;
- int k;
-
- /* Get the absolute distance. */
- ir = 1.0f / sqrtf(r2);
- r = r2 * ir;
-
- /* Evaluate the gravity kernel. */
- kernel_grav_eval(r, &acc);
-
- /* Scale the acceleration. */
- acc *= const_G * ir * ir * ir;
-
- /* Aggregate the accelerations. */
- for (k = 0; k < 3; k++) {
- w = acc * dx[k];
- pi->a_grav[k] -= w * mj;
- pj->a_grav[k] += w * mi;
- }
-}
+ float r2, float *dx, struct gpart *pi, struct gpart *pj) {}
/**
* @brief Gravity potential (Vectorized version)
*/
__attribute__((always_inline)) INLINE static void runner_iact_vec_grav(
- float *R2, float *Dx, struct gpart **pi, struct gpart **pj) {
-
-#ifdef VECTORIZE
-
- vector ir, r, r2, dx[3];
- vector w, acc, ai, aj;
- vector mi, mj;
- int j, k;
-
-#if VEC_SIZE == 8
- mi.v = vec_set(pi[0]->mass, pi[1]->mass, pi[2]->mass, pi[3]->mass,
- pi[4]->mass, pi[5]->mass, pi[6]->mass, pi[7]->mass);
- mj.v = vec_set(pj[0]->mass, pj[1]->mass, pj[2]->mass, pj[3]->mass,
- pj[4]->mass, pj[5]->mass, pj[6]->mass, pj[7]->mass);
- for (k = 0; k < 3; k++)
- dx[k].v = vec_set(Dx[0 + k], Dx[3 + k], Dx[6 + k], Dx[9 + k], Dx[12 + k],
- Dx[15 + k], Dx[18 + k], Dx[21 + k]);
-#elif VEC_SIZE == 4
- mi.v = vec_set(pi[0]->mass, pi[1]->mass, pi[2]->mass, pi[3]->mass);
- mj.v = vec_set(pj[0]->mass, pj[1]->mass, pj[2]->mass, pj[3]->mass);
- for (k = 0; k < 3; k++)
- dx[k].v = vec_set(Dx[0 + k], Dx[3 + k], Dx[6 + k], Dx[9 + k]);
-#endif
-
- /* Get the radius and inverse radius. */
- r2.v = vec_load(R2);
- ir.v = vec_rsqrt(r2.v);
- ir.v = ir.v - vec_set1(0.5f) * ir.v * (r2.v * ir.v * ir.v - vec_set1(1.0f));
- r.v = r2.v * ir.v;
-
- /* Evaluate the gravity kernel. */
- blender_eval_vec(&r, &acc);
-
- /* Scale the acceleration. */
- acc.v *= vec_set1(const_G) * ir.v * ir.v * ir.v;
-
- /* Aggregate the accelerations. */
- for (k = 0; k < 3; k++) {
- w.v = acc.v * dx[k].v;
- ai.v = w.v * mj.v;
- aj.v = w.v * mi.v;
- for (j = 0; j < VEC_SIZE; j++) {
- pi[j]->a_grav[k] -= ai.f[j];
- pj[j]->a_grav[k] += aj.f[j];
- }
- }
-
-#else
-
- for (int k = 0; k < VEC_SIZE; k++)
- runner_iact_grav(R2[k], &Dx[3 * k], pi[k], pj[k]);
-
-#endif
-}
+ float *R2, float *Dx, struct gpart **pi, struct gpart **pj) {}
#endif /* SWIFT_RUNNER_IACT_GRAV_H */
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 2d3bc69ba6f9414a8e9ccaf1b5662227184f7aa1..6eadcd5036d768a413ee475e88f6df4b6e11d64d 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -37,18 +37,17 @@ struct gpart {
/* Particle mass. */
float mass;
+ /* Softening length */
+ float epsilon;
+
/* Particle time of beginning of time-step. */
int ti_begin;
/* Particle time of end of time-step. */
int ti_end;
- /* /\* current time of x, and of v_full *\/ */
- /* float tx; */
- /* float tv; */
-
double mass_interacted;
-
+
/* Anonymous union for id/part. */
union {
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 41f3348db10172c1d0d5813debb50365252a413b..eccfef645f79eb1942a37612cf3490dfff43fb7f 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -83,8 +83,8 @@ __attribute__((always_inline)) INLINE static void hydro_write_particles(
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "SmoothingLength",
FLOAT, N, 1, parts, N_total, mpi_rank, offset, h, us,
UNIT_CONV_LENGTH);
- writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Entropy",
- FLOAT, N, 1, parts, N_total, mpi_rank, offset, entropy, us,
+ writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Entropy", FLOAT, N,
+ 1, parts, N_total, mpi_rank, offset, entropy, us,
UNIT_CONV_ENTROPY_PER_UNIT_MASS);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "ParticleIDs",
ULONGLONG, N, 1, parts, N_total, mpi_rank, offset, id, us,
diff --git a/src/kernel_gravity.c b/src/kernel_gravity.c
index 639a964c813ef7fd95008857ee17b7dd5ffafb27..4a5df7a8772618cf2f5359272a5a0a9f29809e85 100644
--- a/src/kernel_gravity.c
+++ b/src/kernel_gravity.c
@@ -42,7 +42,7 @@ float gadget(float r) {
(21.333333333333 - 48.0 * u + 38.4 * u * u -
10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
}
- return const_G * fac;
+ return fac;
}
/**
@@ -66,7 +66,7 @@ void gravity_kernel_dump(float r_max, int N) {
x4[1] = x4[0];
x4[0] = x;
kernel_grav_eval(x, &w);
- w *= const_G / (x * x * x);
+ w *= 1.f / (x * x * x);
// blender_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
printf(" %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n", x, w * x, w4[0],
w4[1], w4[2], w4[3], gadget(x) * x);
diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h
index cca449efa6dde73eac428d1e20fa5c28156360fa..e976c6108edfb5bcd17dd161b8bb058b03326055 100644
--- a/src/kernel_gravity.h
+++ b/src/kernel_gravity.h
@@ -25,6 +25,7 @@
#include "inline.h"
#include "vector.h"
+#define const_epsilon 1.
#define const_iepsilon (1. / const_epsilon)
#define const_iepsilon2 (const_iepsilon *const_iepsilon)
#define const_iepsilon3 (const_iepsilon2 *const_iepsilon)
diff --git a/src/runner.c b/src/runner.c
index 49223b46a11c1b113b9f0633ccc961e70b9bee43..76f512d9dd5826bf99777829cef0d96a897ece2a 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -725,6 +725,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
const struct external_potential *potential = r->e->external_potential;
const struct hydro_props *hydro_properties = r->e->hydro_properties;
const struct phys_const *constants = r->e->physical_constants;
+ const double const_G = constants->const_newton_G;
const float ln_max_h_change = hydro_properties->log_max_h_change;
const int is_fixdt =
(r->e->policy & engine_policy_fixdt) == engine_policy_fixdt;
@@ -756,7 +757,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
if (is_fixdt || gp->ti_end <= ti_current) {
/* First, finish the force calculation */
- gravity_end_force(gp);
+ gravity_end_force(gp, const_G);
if (gp->id == -ICHECK)
message("id=%lld a=[%f %f %f] M=%f\n", gp->id, gp->a_grav[0],
@@ -853,7 +854,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
/* And do the same of the extra variable */
hydro_end_force(p);
- if (p->gpart != NULL) gravity_end_force(p->gpart);
+ if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
/* Now we are ready to compute the next time-step size */
int new_dti;
diff --git a/src/tools.c b/src/tools.c
index f298b3494da0c77dd76fdfb30940afa3db1a3d8e..a178a2d7d7d7ae5175f1e49f73de6fdcc4ed4a4b 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -490,7 +490,7 @@ void shuffle_particles(struct part *parts, const int count) {
* @brief gparts The array of particles.
* @brief gcount The number of particles.
*/
-void gravity_n2(struct gpart *gparts, const int gcount) {
+void gravity_n2(struct gpart *gparts, const int gcount, const struct phys_const *constants) {
/* Reset everything */
for (int pid = 0; pid < gcount; pid++) {
@@ -533,4 +533,13 @@ void gravity_n2(struct gpart *gparts, const int gcount) {
gpj->a_grav[2] += wdx[2] * mi;
}
}
+
+ /* Multiply by Newton's constant */
+ const double const_G = constants->const_newton_G;
+ for (int pid = 0; pid < gcount; pid++) {
+ struct gpart *restrict gpi = &gparts[pid];
+ gpi->a_grav[0] *= const_G;
+ gpi->a_grav[1] *= const_G;
+ gpi->a_grav[2] *= const_G;
+ }
}
diff --git a/src/tools.h b/src/tools.h
index 90bc4da01134b9f929ecb5a0f36538be30e995e1..4f587b2ab1211ebe2c329e0b7362b7ade005d2a7 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -22,8 +22,9 @@
#ifndef SWIFT_TOOL_H
#define SWIFT_TOOL_H
-#include "runner.h"
#include "cell.h"
+#include "physical_constants.h"
+#include "runner.h"
void factor(int value, int *f1, int *f2);
void density_dump(int N);
@@ -40,6 +41,7 @@ void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
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);
+void gravity_n2(struct gpart *gparts, const int gcount,
+ const struct phys_const *constants);
#endif /* SWIFT_TOOL_H */