diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index ce9e461fdce19bb57c071d577f3271fce000e916..fad0700cf18f5650b393e41028e54d69ca305137 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -20,19 +20,25 @@
#include <float.h>
#include "potentials.h"
+
/**
* @brief Computes the gravity time-step of a given particle.
*
+ * This function only branches towards the potential chosen by the user.
+ *
+ * @param potential The properties of the external potential.
* @param phys_cont The physical constants in internal units.
* @param gp Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static float gravity_compute_timestep(
- const struct external_potential* potential, const struct phys_const* const phys_const, const struct gpart* const gp) {
+ const struct external_potential* potential,
+ const struct phys_const* const phys_const, const struct gpart* const gp) {
float dt = FLT_MAX;
#ifdef EXTERNAL_POTENTIAL_POINTMASS
- dt = fminf(dt, external_gravity_pointmass_timestep(potential, phys_const, gp));
+ dt =
+ fminf(dt, external_gravity_pointmass_timestep(potential, phys_const, gp));
#endif
return dt;
@@ -79,10 +85,15 @@ __attribute__((always_inline))
/**
* @brief Computes the gravitational acceleration induced by external potentials
*
+ * This function only branches towards the potential chosen by the user.
+ *
+ * @param potential The properties of the external potential.
* @param phys_const The physical constants in internal units.
* @param gp The particle to act upon.
*/
-__attribute__((always_inline)) INLINE static void external_gravity(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* gp) {
+__attribute__((always_inline)) INLINE static void external_gravity(
+ const struct external_potential* potential,
+ const struct phys_const* const phys_const, struct gpart* gp) {
#ifdef EXTERNAL_POTENTIAL_POINTMASS
external_gravity_pointmass(potential, phys_const, gp);
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
index 654745bfeb70dddba772af9e23797713376377a7..72ae5fd55f399448d78c6f2d1533c7f007f16dd1 100644
--- a/src/gravity/Default/gravity_debug.h
+++ b/src/gravity/Default/gravity_debug.h
@@ -24,5 +24,6 @@ __attribute__((always_inline))
"v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"mass=%.3e t_begin=%d, t_end=%d\n",
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], p->mass, p->ti_begin, p->ti_end);
+ p->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->ti_begin,
+ p->ti_end);
}
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 6a33f786cd2ee1bb6fdbba20e544fa4b9d2dd556..7f7990ef11f6fe65b39b011e4a646b595911c862 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -40,10 +40,9 @@ struct gpart {
/* Particle time of end of time-step. */
int ti_end;
-
- /* current time of x, and of v_full */
- float tx;
- float tv;
+ /* /\* current time of x, and of v_full *\/ */
+ /* float tx; */
+ /* float tv; */
/* Anonymous union for id/part. */
union {
diff --git a/src/runner.c b/src/runner.c
index 331b91d4380fe6925f1cad527f3fb16c9698716e..512fdbdc180f2764b3def4a7bfb379d45360f413 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -71,19 +71,20 @@ const float runner_shift[13 * 3] = {
/* Does the axis need flipping ? */
const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
-#define MY_CELL 428428428
-#define DX 0.1
-#define NX 1000
-#define CELL_ID \
- ((int)(c->loc[0] / DX) * NX *NX + (int)(c->loc[1] / DX) * NX + \
- (int)(c->loc[2] / DX))
-#define OUT \
- if (CELL_ID == MY_CELL) { \
- message(" cell= %d gcount=%d time=%f \n", CELL_ID, c->gcount, r->e->time); \
- }
-//#define OUT message(" cell %d %d %f \n",CELL_ID, c->count, r->e->time);
-//#define OUT if(CELL_ID == MY_CELL) message("\n cell %f %f %f %d %d
-//%f\n",c->loc[0],c->loc[1],c->loc[2], CELL_ID, c->count, r->e->time);
+
+/* #define MY_CELL 428428428 */
+/* #define DX 0.1 */
+/* #define NX 1000 */
+/* #define CELL_ID \ */
+/* ((int)(c->loc[0] / DX) * NX *NX + (int)(c->loc[1] / DX) * NX + \ */
+/* (int)(c->loc[2] / DX)) */
+/* #define OUT \ */
+/* if (CELL_ID == MY_CELL) { \ */
+/* message(" cell= %d gcount=%d time=%f \n", CELL_ID, c->gcount, r->e->time); \ */
+/* } */
+/* #define OUT message(" cell %d %d %f \n",CELL_ID, c->count, r->e->time); */
+/* #define OUT if(CELL_ID == MY_CELL) message("\n cell %f %f %f %d %d */
+/* %f\n",c->loc[0],c->loc[1],c->loc[2], CELL_ID, c->count, r->e->time); */
/* Import the density loop functions. */
#define FUNCTION density
@@ -106,20 +107,11 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
void runner_dograv_external(struct runner *r, struct cell *c) {
struct gpart *g, *gparts = c->gparts;
- float L[3], E;
int i, k, gcount = c->gcount;
const int ti_current = r->e->ti_current;
- // struct space *s = r->e->s;
- // double CentreOfPotential[3];
TIMER_TIC;
- /* /\* location of external gravity point mass - should pass in as
- * paraneter *\/ */
- /* CentreOfPotential[0] = 0.5 * s->dim[0]; */
- /* CentreOfPotential[1] = 0.5 * s->dim[1]; */
- /* CentreOfPotential[2] = 0.5 * s->dim[2]; */
-
/* Recurse? */
if (c->split) {
for (k = 0; k < 8; k++)
@@ -130,6 +122,7 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
#ifdef TASK_VERBOSE
OUT;
#endif
+
/* Loop over the parts in this cell. */
for (i = 0; i < gcount; i++) {
@@ -138,49 +131,50 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
/* Is this part within the time step? */
if (g->ti_end <= ti_current) {
- // external_gravity_pointmass(e->physical_constants,
- // potential_constants, g);
+
external_gravity(r->e->potential, r->e->physical_constants, g);
- /* check for energy and angular momentum conservation - begin by
- * synchronizing velocity*/
- const float dx = g->x[0] - r->e->potential->point_mass.x;
- const float dy = g->x[1] - r->e->potential->point_mass.y;
- const float dz = g->x[2] - r->e->potential->point_mass.z;
- const float dr = sqrtf((dx * dx) + (dy * dy) + (dz * dz));
- const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
-
- const int current_dti = g->ti_end - g->ti_begin;
- const float dt = 0.5f * current_dti * r->e->timeBase;
- const float vx = g->v_full[0] + dt * g->a_grav[0];
- const float vy = g->v_full[1] + dt * g->a_grav[1];
- const float vz = g->v_full[2] + dt * g->a_grav[2];
-
- /* E/L */
- E = 0.5 * ((vx * vx) + (vy * vy) + (vz * vz)) -
- r->e->physical_constants->newton_gravity *
- r->e->potential->point_mass.mass * rinv;
- L[0] = dy * vz - dz * vy;
- L[1] = dz * vx - dx * vz;
- L[2] = dx * vy - dy * vx;
- if (abs(g->id) == 1) {
- float v2 = vx * vx + vy * vy + vz * vz;
- float fg = r->e->physical_constants->newton_gravity *
- r->e->potential->point_mass.mass * rinv;
- float fga = sqrtf((g->a_grav[0] * g->a_grav[0]) +
- (g->a_grav[1] * g->a_grav[1]) +
- (g->a_grav[2] * g->a_grav[2])) *
- dr;
- // message("grav_external time= %f\t V_c^2= %f GM/r= %f E= %f L[2]= %f
- // x= %f y= %f vx= %f vy= %f\n", r->e->time, v2, fg, E, L[2], g->x[0],
- // g->x[1], vx, vy);
- message("%f\t %f %f %f %f %f %f %f %f %f %f %f %f %f\n", r->e->time,
- g->tx, g->tv, dt, v2, fg, fga, dr, E, L[2], g->x[0], g->x[1],
- vx, vy);
- /* message(" G=%e M=%e\n", r->e->physical_constants->newton_gravity,
- * r->e->potential->point_mass.mass); */
- /* exit(-1); */
- }
+ /* /\* check for energy and angular momentum conservation - begin by */
+ /* * synchronizing velocity*\/ */
+ /* const float dx = g->x[0] - r->e->potential->point_mass.x; */
+ /* const float dy = g->x[1] - r->e->potential->point_mass.y; */
+ /* const float dz = g->x[2] - r->e->potential->point_mass.z; */
+ /* const float dr = sqrtf((dx * dx) + (dy * dy) + (dz * dz)); */
+ /* const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz); */
+
+ /* const int current_dti = g->ti_end - g->ti_begin; */
+ /* const float dt = 0.5f * current_dti * r->e->timeBase; */
+ /* const float vx = g->v_full[0] + dt * g->a_grav[0]; */
+ /* const float vy = g->v_full[1] + dt * g->a_grav[1]; */
+ /* const float vz = g->v_full[2] + dt * g->a_grav[2]; */
+
+ /* /\* E/L *\/ */
+ /* float L[3], E; */
+ /* E = 0.5 * ((vx * vx) + (vy * vy) + (vz * vz)) - */
+ /* r->e->physical_constants->newton_gravity * */
+ /* r->e->potential->point_mass.mass * rinv; */
+ /* L[0] = dy * vz - dz * vy; */
+ /* L[1] = dz * vx - dx * vz; */
+ /* L[2] = dx * vy - dy * vx; */
+ /* if (abs(g->id) == 1) { */
+ /* float v2 = vx * vx + vy * vy + vz * vz; */
+ /* float fg = r->e->physical_constants->newton_gravity * */
+ /* r->e->potential->point_mass.mass * rinv; */
+ /* float fga = sqrtf((g->a_grav[0] * g->a_grav[0]) + */
+ /* (g->a_grav[1] * g->a_grav[1]) + */
+ /* (g->a_grav[2] * g->a_grav[2])) * */
+ /* dr; */
+ /* // message("grav_external time= %f\t V_c^2= %f GM/r= %f E= %f L[2]= %f */
+ /* // x= %f y= %f vx= %f vy= %f\n", r->e->time, v2, fg, E, L[2], g->x[0], */
+ /* // g->x[1], vx, vy); */
+ /* message("%f\t %f %f %f %f %f %f %f %f %f %f %f %f %f\n", r->e->time, */
+ /* g->tx, g->tv, dt, v2, fg, fga, dr, E, L[2], g->x[0], g->x[1], */
+ /* vx, vy); */
+ /* /\* message(" G=%e M=%e\n", r->e->physical_constants->newton_gravity, */
+ /* * r->e->potential->point_mass.mass); *\/ */
+ /* /\* exit(-1); *\/ */
+ /* } */
+
}
}
TIMER_TOC(timer_dograv_external);
diff --git a/src/task.h b/src/task.h
index 56dff013356884543948114d6e5d07c4fb434c4e..25cc886f4b38456a0431fb6c7d0b7b1864053dd9 100644
--- a/src/task.h
+++ b/src/task.h
@@ -28,7 +28,6 @@
#include "cycle.h"
/* Some constants. */
-#define TASK_VERBOSE
#define task_maxwait 3
#define task_maxunlock 15