diff --git a/examples/main.c b/examples/main.c
index 2059a8018ba5b51e65ef78a1c1e40027a55af671..38db83a334f6cbdc04afed45a12446dd82e84f1d 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -88,8 +88,8 @@ void print_help_message() {
int main(int argc, char *argv[]) {
struct clocks_time tic, toc;
- int nr_nodes = 1, myrank = 0;
+ int nr_nodes = 1, myrank = 0;
#ifdef WITH_MPI
/* Start by initializing MPI. */
int res = 0, prov = 0;
@@ -258,14 +258,28 @@ int main(int argc, char *argv[]) {
MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
#endif
- /* Initialize unit system */
+/* Prepare the domain decomposition scheme */
+#ifdef WITH_MPI
+ struct partition initial_partition;
+ enum repartition_type reparttype;
+ partition_init(&initial_partition, &reparttype, params, nr_nodes);
+
+ /* Let's report what we did */
+ if (myrank == 0) {
+ message("Using initial partition %s",
+ initial_partition_name[initial_partition.type]);
+ if (initial_partition.type == INITPART_GRID)
+ message("grid set to [ %i %i %i ].", initial_partition.grid[0],
+ initial_partition.grid[1], initial_partition.grid[2]);
+ message("Using %s repartitioning", repartition_name[reparttype]);
+ }
+#endif
+
+ /* Initialize unit system and constants */
struct UnitSystem us;
struct phys_const prog_const;
- struct external_potential potential;
units_init(&us, params);
initPhysicalConstants(&us, &prog_const);
- if( with_external_gravity)
- initPotentialProperties(params, &us, &potential);
if (myrank == 0) {
message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
message("Unit system: U_L = %e cm.", us.UnitLength_in_cgs);
@@ -274,31 +288,18 @@ int main(int argc, char *argv[]) {
message("Unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
message("Density units: %e a^%f h^%f.",
units_conversion_factor(&us, UNIT_CONV_DENSITY),
- units_a_factor(&us, UNIT_CONV_DENSITY),
+ units_a_factor(&us, UNIT_CONV_DENSITY),
units_h_factor(&us, UNIT_CONV_DENSITY));
message("Entropy units: %e a^%f h^%f.",
units_conversion_factor(&us, UNIT_CONV_ENTROPY),
- units_a_factor(&us, UNIT_CONV_ENTROPY),
+ units_a_factor(&us, UNIT_CONV_ENTROPY),
units_h_factor(&us, UNIT_CONV_ENTROPY));
message("Gravity constant = %e", prog_const.newton_gravity);
}
-/* Prepare the domain decomposition scheme */
-#ifdef WITH_MPI
- struct partition initial_partition;
- enum repartition_type reparttype;
- partition_init(&initial_partition, &reparttype, params, nr_nodes);
-
- /* Let's report what we did */
- if (myrank == 0) {
- message("Using initial partition %s",
- initial_partition_name[initial_partition.type]);
- if (initial_partition.type == INITPART_GRID)
- message("grid set to [ %i %i %i ].", initial_partition.grid[0],
- initial_partition.grid[1], initial_partition.grid[2]);
- message("Using %s repartitioning", repartition_name[reparttype]);
- }
-#endif
+ /* Initialise the external potential properties */
+ struct external_potential potential;
+ if (with_external_gravity) initPotentialProperties(params, &us, &potential);
/* Read particles and space information from (GADGET) ICs */
char ICfileName[200] = "";
@@ -475,16 +476,15 @@ int main(int argc, char *argv[]) {
int icount = 0;
space_map_cells_pre(&s, 0, &map_cellcheck, &icount);
message("map_cellcheck picked up %i parts.", icount);
-
}
if (myrank == 0) {
- int data[2] = {s.maxdepth ,0};
+ int data[2] = {s.maxdepth, 0};
space_map_cells_pre(&s, 0, &map_maxdepth, data);
message("nr of cells at depth %i is %i.", data[0], data[1]);
}
- //exit(-99);
+ // exit(-99);
/* Legend */
if (myrank == 0)
printf("# %6s %14s %14s %10s %10s %16s [%s]\n", "Step", "Time", "Time-step",
diff --git a/src/engine.c b/src/engine.c
index b9d9f4c48f9e0e1b4027ba240657a25e7ff8e439..364370f73cf690f26d64ce1e1331122190645bb5 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1895,7 +1895,7 @@ void engine_init_particles(struct engine *e) {
/* We always have sort tasks */
mask |= 1 << task_type_sort;
mask |= 1 << task_type_init;
-
+
/* Add the tasks corresponding to hydro operations to the masks */
if ((e->policy & engine_policy_hydro) == engine_policy_hydro) {
@@ -2343,7 +2343,7 @@ static bool hyperthreads_present(void) {
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
- int policy, int verbose,
+ int policy, int verbose,
const struct phys_const *physical_constants,
const struct external_potential *potential) {
diff --git a/src/physical_constants_cgs.h b/src/physical_constants_cgs.h
index 575ebffb7897a90c5b3a016a4d800a3c827d33f9..30bee1367a26adff8f1bc491cc5754d103053db8 100644
--- a/src/physical_constants_cgs.h
+++ b/src/physical_constants_cgs.h
@@ -23,7 +23,7 @@
/* physical constants in cgs */
#define NEWTON_GRAVITY_CGS 6.672e-8
#define SOLAR_MASS_IN_CGS 1.989e33
-#define PARSEC_IN_CGS 3.0856776e18
+#define PARSEC_IN_CGS 3.0856776e18
#define PROTON_MASS_IN_CGS 1.6726231e24
#define YEAR_IN_CGS 3.154e+7
diff --git a/src/potentials.c b/src/potentials.c
index 324dfcad74f5ecb74fe0ed6ba3c4ffbac81917eb..6cb2ab6188944f0889bdbfbb4a314b2084c70b06 100644
--- a/src/potentials.c
+++ b/src/potentials.c
@@ -31,21 +31,31 @@
* @param us The current internal system of units
* @param potential The external potential properties to initialize
*/
-void initPotentialProperties(const struct swift_params * parameter_file,
- struct UnitSystem* us,
+void initPotentialProperties(const struct swift_params* parameter_file,
+ struct UnitSystem* us,
struct external_potential* potential) {
- potential->point_mass.x = parser_get_param_double(parameter_file, "PointMass:position_x");
- potential->point_mass.y = parser_get_param_double(parameter_file, "PointMass:position_y");
- potential->point_mass.z = parser_get_param_double(parameter_file, "PointMass:position_z");
- potential->point_mass.mass = parser_get_param_double(parameter_file, "PointMass:mass");
- message("Point mass properties are (x,y,z) = (%e, %e, %e), M = %e", potential->point_mass.x, potential->point_mass.y, potential->point_mass.z, potential->point_mass.mass);
-
+ potential->point_mass.x =
+ parser_get_param_double(parameter_file, "PointMass:position_x");
+ potential->point_mass.y =
+ parser_get_param_double(parameter_file, "PointMass:position_y");
+ potential->point_mass.z =
+ parser_get_param_double(parameter_file, "PointMass:position_z");
+ potential->point_mass.mass =
+ parser_get_param_double(parameter_file, "PointMass:mass");
+ message("Point mass properties are (x,y,z) = (%e, %e, %e), M = %e",
+ potential->point_mass.x, potential->point_mass.y,
+ potential->point_mass.z, potential->point_mass.mass);
+
/* potential->point_mass.x = */
- /* 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH); */
+ /* 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH);
+ */
/* potential->point_mass.y = */
- /* 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH); */
+ /* 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH);
+ */
/* potential->point_mass.z = */
- /* 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH); */
+ /* 50000 * PARSEC_IN_CGS / units_conversion_factor(us, UNIT_CONV_LENGTH);
+ */
/* potential->point_mass.mass = */
- /* 1e10 * SOLAR_MASS_IN_CGS / units_conversion_factor(us, UNIT_CONV_MASS); */
+ /* 1e10 * SOLAR_MASS_IN_CGS / units_conversion_factor(us, UNIT_CONV_MASS);
+ */
}
diff --git a/src/potentials.h b/src/potentials.h
index e31ac207974669347541439479922d2347be163b..e489781a62890bf880aee8003df3109938278c11 100644
--- a/src/potentials.h
+++ b/src/potentials.h
@@ -42,9 +42,12 @@ struct external_potential {
#ifdef EXTERNAL_POTENTIAL_POINTMASS
/* #define const_unit_length_in_cgs (1e3 * PARSEC_IN_CGS) */
/* #define const_unit_mass_in_cgs (SOLAR_MASS_IN_CGS) */
-/* #define External_Potential_X (50000 * PARSEC_IN_CGS / const_unit_length_in_cgs) */
-/* #define External_Potential_Y (50000 * PARSEC_IN_CGS / const_unit_length_in_cgs) */
-/* #define External_Potential_Z (50000 * PARSEC_IN_CGS / const_unit_length_in_cgs) */
+/* #define External_Potential_X (50000 * PARSEC_IN_CGS /
+ * const_unit_length_in_cgs) */
+/* #define External_Potential_Y (50000 * PARSEC_IN_CGS /
+ * const_unit_length_in_cgs) */
+/* #define External_Potential_Z (50000 * PARSEC_IN_CGS /
+ * const_unit_length_in_cgs) */
/* #define External_Potential_Mass \ */
/* (1e10 * SOLAR_MASS_IN_CGS / const_unit_mass_in_cgs) */
@@ -56,7 +59,7 @@ struct external_potential {
*/
__attribute__((always_inline))
INLINE static float external_gravity_pointmass_timestep(
- const struct external_potential* potential,
+ const struct external_potential* potential,
const struct phys_const* const phys_const,
const struct gpart* const g) {
@@ -68,12 +71,12 @@ __attribute__((always_inline))
const float drdv = (g->x[0] - potential->point_mass.x) * (g->v_full[0]) +
(g->x[1] - potential->point_mass.y) * (g->v_full[1]) +
(g->x[2] - potential->point_mass.z) * (g->v_full[2]);
- const float dota_x = G_newton * potential->point_mass.mass * rinv * rinv * rinv *
- (-g->v_full[0] + 3.f * rinv * rinv * drdv * dx);
- const float dota_y = G_newton * potential->point_mass.mass * rinv * rinv * rinv *
- (-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
- const float dota_z = G_newton * potential->point_mass.mass * rinv * rinv * rinv *
- (-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
+ const float dota_x = G_newton * potential->point_mass.mass * rinv * rinv *
+ rinv * (-g->v_full[0] + 3.f * rinv * rinv * drdv * dx);
+ const float dota_y = G_newton * potential->point_mass.mass * rinv * rinv *
+ rinv * (-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
+ const float dota_z = G_newton * potential->point_mass.mass * rinv * rinv *
+ rinv * (-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
g->a_grav[2] * g->a_grav[2];
@@ -88,8 +91,9 @@ __attribute__((always_inline))
* @param phys_cont The physical constants in internal units.
* @param gp Pointer to the g-particle data.
*/
-__attribute__((always_inline)) INLINE static void external_gravity_pointmass(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
-
+__attribute__((always_inline)) INLINE static void external_gravity_pointmass(
+ const struct external_potential* potential,
+ const struct phys_const* const phys_const, struct gpart* g) {
/* message(" (x,y,z) = (%e, %e, %e), M= %e", potential->point_mass.x, */
/* potential->point_mass.y, potential->point_mass.z, */
@@ -102,9 +106,12 @@ __attribute__((always_inline)) INLINE static void external_gravity_pointmass(con
const float dz = g->x[2] - potential->point_mass.z;
const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
- g->a_grav[0] += -G_newton * potential->point_mass.mass * dx * rinv * rinv * rinv;
- g->a_grav[1] += -G_newton * potential->point_mass.mass * dy * rinv * rinv * rinv;
- g->a_grav[2] += -G_newton * potential->point_mass.mass * dz * rinv * rinv * rinv;
+ g->a_grav[0] +=
+ -G_newton * potential->point_mass.mass * dx * rinv * rinv * rinv;
+ g->a_grav[1] +=
+ -G_newton * potential->point_mass.mass * dy * rinv * rinv * rinv;
+ g->a_grav[2] +=
+ -G_newton * potential->point_mass.mass * dz * rinv * rinv * rinv;
}
#endif /* EXTERNAL_POTENTIAL_POINTMASS */
@@ -115,8 +122,8 @@ __attribute__((always_inline)) INLINE static void external_gravity_pointmass(con
* @param us The current internal system of units
* @param potential The external potential properties to initialize
*/
-void initPotentialProperties(const struct swift_params * parameter_file,
- struct UnitSystem* us,
+void initPotentialProperties(const struct swift_params* parameter_file,
+ struct UnitSystem* us,
struct external_potential* potential);
#endif /* SWIFT_POTENTIALS_H */
diff --git a/src/runner.c b/src/runner.c
index aab2ad01697e6655c45768ceb6e45cad66229888..331b91d4380fe6925f1cad527f3fb16c9698716e 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -120,7 +120,6 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
/* CentreOfPotential[1] = 0.5 * s->dim[1]; */
/* CentreOfPotential[2] = 0.5 * s->dim[2]; */
-
/* Recurse? */
if (c->split) {
for (k = 0; k < 8; k++)
@@ -151,7 +150,6 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
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];
@@ -160,13 +158,15 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
/* E/L */
E = 0.5 * ((vx * vx) + (vy * vy) + (vz * vz)) -
- r->e->physical_constants->newton_gravity * r->e->potential->point_mass.mass * rinv;
+ 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 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])) *
@@ -174,10 +174,12 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
// 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); */
+ 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); */
}
}
}
@@ -940,7 +942,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
} else {
/* Compute the next timestep (gravity condition) */
- float new_dt = gravity_compute_timestep(r->e->potential, r->e->physical_constants, gp);
+ float new_dt = gravity_compute_timestep(r->e->potential,
+ r->e->physical_constants, gp);
/* Limit timestep within the allowed range */
new_dt = fminf(new_dt, global_dt_max);
@@ -1025,8 +1028,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* Compute the next timestep (gravity condition) */
float new_dt_grav = FLT_MAX;
if (p->gpart != NULL)
- new_dt_grav =
- gravity_compute_timestep(r->e->potential, r->e->physical_constants, p->gpart);
+ new_dt_grav = gravity_compute_timestep(
+ r->e->potential, r->e->physical_constants, p->gpart);
float new_dt = fminf(new_dt_hydro, new_dt_grav);
diff --git a/src/space.c b/src/space.c
index 66bc165b4a2a72a1f553022b6d8905999b122a0b..2e4fe71bfd248948ffd54eb17782269bfba178fb 100644
--- a/src/space.c
+++ b/src/space.c
@@ -1348,9 +1348,9 @@ void space_init(struct space *s, const struct swift_params *params,
space_maxsize = parser_get_param_int(params, "Scheduler:cell_max_size");
space_subsize = parser_get_param_int(params, "Scheduler:cell_sub_size");
space_splitsize = parser_get_param_int(params, "Scheduler:cell_split_size");
- if(verbose)
+ if (verbose)
message("max_size set to %d, sub_size set to %d, split_size set to %d",
- space_maxsize, space_subsize, space_splitsize);
+ space_maxsize, space_subsize, space_splitsize);
/* Check that we have enough cells */
if (s->cell_min * 3 > dim[0] || s->cell_min * 3 > dim[1] ||
diff --git a/src/task.c b/src/task.c
index ccf408ad770ba8cdd46a7f1d8bbc762a51b9c523..31c2bfcb3462b9f91547ebd631f645e190c07143 100644
--- a/src/task.c
+++ b/src/task.c
@@ -47,9 +47,9 @@
/* Task type names. */
const char *taskID_names[task_type_count] = {
- "none", "sort", "self", "pair", "sub",
- "init", "ghost", "drift", "kick", "send",
- "recv", "grav_pp", "grav_mm", "grav_up", "grav_down",
+ "none", "sort", "self", "pair", "sub",
+ "init", "ghost", "drift", "kick", "send",
+ "recv", "grav_pp", "grav_mm", "grav_up", "grav_down",
"grav_external", "part_sort", "gpart_sort", "split_cell", "rewait"};
const char *subtaskID_names[task_type_count] = {"none", "density",