diff --git a/examples/ExternalPointMass/run.sh b/examples/ExternalPointMass/run.sh
index 9f90ca395a5c8cf83e67928b3fdbd4d8529ac254..4ac513f09cb8ac8dcefc256a68478e215b8bc320 100755
--- a/examples/ExternalPointMass/run.sh
+++ b/examples/ExternalPointMass/run.sh
@@ -7,4 +7,4 @@ then
python makeIC.py 10000
fi
-../swift -g -t 2 externalPointMass.yml 2>&1 | tee output.log
+../swift -g -t 1 externalPointMass.yml 2>&1 | tee output.log
diff --git a/examples/IsothermalPotential/GravityOnly/README b/examples/IsothermalPotential/GravityOnly/README
index 90fb1872aa2301cab133b8a20a8cd8de724d4553..1621aaa8ab90420dcf69b5b9caea394d619b62cf 100644
--- a/examples/IsothermalPotential/GravityOnly/README
+++ b/examples/IsothermalPotential/GravityOnly/README
@@ -1,6 +1,3 @@
-;
-; this probelm generates a set of gravity particles in an isothermal
-; potential and follows their orbits. Tests verify consdevation of
-; energy and angular momentum
-;
-;
+This example generates a set of particles in an isothermal potential
+and follows their orbits. IDL scripts verify the conservation of
+energy and angular momentum.
diff --git a/src/cell.c b/src/cell.c
index a0da9bf58ddad6fad48b7146a4b97fd7c800cc50..5d3f0673f77aca351459177ebecddda54106c6a6 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1007,3 +1007,23 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
return 0;
}
+
+/**
+ * @brief Set the super-cell pointers for all cells in a hierarchy.
+ *
+ * @param c The top-level #cell to play with.
+ * @param super Pointer to the deepest cell with tasks in this part of the tree.
+ */
+void cell_set_super(struct cell *c, struct cell *super) {
+
+ /* Are we in a cell with some kind of self/pair task ? */
+ if (super == NULL && c->nr_tasks > 0) super = c;
+
+ /* Set the super-cell */
+ c->super = super;
+
+ /* Recurse */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) cell_set_super(c->progeny[k], super);
+}
diff --git a/src/cell.h b/src/cell.h
index 4b970d587923f53dab909250de2b31da3252d6ae..7c6eed55b83277fe9c489d2ab58683d4859d3d59 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -131,14 +131,9 @@ struct cell {
/* Parent cell. */
struct cell *parent;
- /* Super cell, i.e. the highest-level supercell that has hydro interactions.
- */
+ /* Super cell, i.e. the highest-level supercell that has pair/self tasks */
struct cell *super;
- /* Super cell, i.e. the highest-level supercell that has gravity interactions.
- */
- struct cell *gsuper;
-
/* The task computing this cell's sorts. */
struct task *sorts;
int sortsize;
@@ -242,5 +237,6 @@ void cell_check_multipole(struct cell *c, void *data);
void cell_clean(struct cell *c);
int cell_is_drift_needed(struct cell *c, int ti_current);
int cell_unskip_tasks(struct cell *c, struct scheduler *s);
+void cell_set_super(struct cell *c, struct cell *super);
#endif /* SWIFT_CELL_H */
diff --git a/src/engine.c b/src/engine.c
index a41659ada346d44d105b65a678995fab1474ce78..a8773fcbeab48bc775503c1795508f5a9bcca57f 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -111,119 +111,57 @@ void engine_addlink(struct engine *e, struct link **l, struct task *t) {
res->next = atomic_swap(l, res);
}
-/**
- * @brief Generate the gravity hierarchical tasks for a hierarchy of cells -
- * i.e. all the O(Npart) tasks.
- *
- * Tasks are only created here. The dependencies will be added later on.
- *
- * @param e The #engine.
- * @param c The #cell.
- * @param gsuper The gsuper #cell.
- */
-void engine_make_gravity_hierarchical_tasks(struct engine *e, struct cell *c,
- struct cell *gsuper) {
-
- struct scheduler *s = &e->sched;
- const int is_with_external_gravity =
- (e->policy & engine_policy_external_gravity) ==
- engine_policy_external_gravity;
- const int is_fixdt = (e->policy & engine_policy_fixdt) == engine_policy_fixdt;
-
- /* Is this the super-cell? */
- if (gsuper == NULL && (c->grav != NULL || (!c->split && c->gcount > 0))) {
-
- /* This is the super cell, i.e. the first with gravity tasks attached. */
- gsuper = c;
-
- /* Local tasks only... */
- if (c->nodeID == e->nodeID) {
-
- /* Add the init task. */
- if (c->init == NULL)
- c->init = scheduler_addtask(s, task_type_init, task_subtype_none, 0, 0,
- c, NULL, 0);
-
- /* Add the kick task that matches the policy. */
- if (is_fixdt) {
- if (c->kick == NULL)
- c->kick = scheduler_addtask(s, task_type_kick_fixdt,
- task_subtype_none, 0, 0, c, NULL, 0);
- } else {
- if (c->kick == NULL)
- c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0,
- 0, c, NULL, 0);
- }
-
- if (is_with_external_gravity)
- c->grav_external = scheduler_addtask(
- s, task_type_grav_external, task_subtype_none, 0, 0, c, NULL, 0);
- }
- }
-
- /* Set the super-cell. */
- c->gsuper = gsuper;
-
- /* Recurse. */
- if (c->split)
- for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL)
- engine_make_gravity_hierarchical_tasks(e, c->progeny[k], gsuper);
-}
-
/**
* @brief Generate the hydro hierarchical tasks for a hierarchy of cells -
* i.e. all the O(Npart) tasks.
*
* Tasks are only created here. The dependencies will be added later on.
*
+ * Note that there is no need to recurse below the super-cell.
+ *
* @param e The #engine.
* @param c The #cell.
- * @param super The super #cell.
*/
-void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
- struct cell *super) {
+void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
struct scheduler *s = &e->sched;
- const int is_fixdt = (e->policy & engine_policy_fixdt) == engine_policy_fixdt;
- const int is_with_cooling =
- (e->policy & engine_policy_cooling) == engine_policy_cooling;
- const int is_with_sourceterms =
- (e->policy & engine_policy_sourceterms) == engine_policy_sourceterms;
+ const int is_fixdt = (e->policy & engine_policy_fixdt);
+ const int is_hydro = (e->policy & engine_policy_hydro);
+ const int is_with_cooling = (e->policy & engine_policy_cooling);
+ const int is_with_sourceterms = (e->policy & engine_policy_sourceterms);
- /* Is this the super-cell? */
- if (super == NULL && (c->density != NULL || (c->count > 0 && !c->split))) {
-
- /* This is the super cell, i.e. the first with density tasks attached. */
- super = c;
+ /* Are we in a super-cell ? */
+ if (c->super == c) {
/* Local tasks only... */
if (c->nodeID == e->nodeID) {
/* Add the init task. */
- if (c->init == NULL)
- c->init = scheduler_addtask(s, task_type_init, task_subtype_none, 0, 0,
- c, NULL, 0);
+ c->init = scheduler_addtask(s, task_type_init, task_subtype_none, 0, 0, c,
+ NULL, 0);
/* Add the kick task that matches the policy. */
if (is_fixdt) {
- if (c->kick == NULL)
- c->kick = scheduler_addtask(s, task_type_kick_fixdt,
- task_subtype_none, 0, 0, c, NULL, 0);
+ c->kick = scheduler_addtask(s, task_type_kick_fixdt, task_subtype_none,
+ 0, 0, c, NULL, 0);
} else {
- if (c->kick == NULL)
- c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0,
- 0, c, NULL, 0);
+ c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0, 0,
+ c, NULL, 0);
}
/* Generate the ghost task. */
- c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0, 0,
- c, NULL, 0);
+ if (is_hydro)
+ c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0,
+ 0, c, NULL, 0);
+
#ifdef EXTRA_HYDRO_LOOP
/* Generate the extra ghost task. */
- c->extra_ghost = scheduler_addtask(s, task_type_extra_ghost,
- task_subtype_none, 0, 0, c, NULL, 0);
+ if (is_hydro)
+ c->extra_ghost = scheduler_addtask(s, task_type_extra_ghost,
+ task_subtype_none, 0, 0, c, NULL, 0);
#endif
+
+ /* Cooling task */
if (is_with_cooling)
c->cooling = scheduler_addtask(s, task_type_cooling, task_subtype_none,
0, 0, c, NULL, 0);
@@ -232,16 +170,19 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
c->sourceterms = scheduler_addtask(s, task_type_sourceterms,
task_subtype_none, 0, 0, c, NULL, 0);
}
- }
- /* Set the super-cell. */
- c->super = super;
+ } else { /* We are above the super-cell so need to go deeper */
- /* Recurse. */
- if (c->split)
- for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL)
- engine_make_hydro_hierarchical_tasks(e, c->progeny[k], super);
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->super != NULL) error("Incorrectly set super pointer");
+#endif
+
+ /* Recurse. */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ engine_make_hierarchical_tasks(e, c->progeny[k]);
+ }
}
/**
@@ -1289,6 +1230,30 @@ void engine_make_gravity_tasks(struct engine *e) {
}
}
+void engine_make_external_gravity_tasks(struct engine *e) {
+
+ struct space *s = e->s;
+ struct scheduler *sched = &e->sched;
+ const int nodeID = e->nodeID;
+ struct cell *cells = s->cells_top;
+ const int nr_cells = s->nr_cells;
+
+ for (int cid = 0; cid < nr_cells; ++cid) {
+
+ struct cell *ci = &cells[cid];
+
+ /* Skip cells without gravity particles */
+ if (ci->gcount == 0) continue;
+
+ /* Is that neighbour local ? */
+ if (ci->nodeID != nodeID) continue;
+
+ /* If the cells is local build a self-interaction */
+ scheduler_addtask(sched, task_type_self, task_subtype_external_grav, 0, 0,
+ ci, NULL, 0);
+ }
+}
+
/**
* @brief Constructs the top-level pair tasks for the first hydro loop over
* neighbours
@@ -1431,11 +1396,27 @@ static inline void engine_make_gravity_dependencies(struct scheduler *sched,
struct cell *c) {
/* init --> gravity --> kick */
- scheduler_addunlock(sched, c->gsuper->init, gravity);
- scheduler_addunlock(sched, gravity, c->gsuper->kick);
+ scheduler_addunlock(sched, c->super->init, gravity);
+ scheduler_addunlock(sched, gravity, c->super->kick);
/* grav_up --> gravity ( --> kick) */
- scheduler_addunlock(sched, c->gsuper->grav_up, gravity);
+ scheduler_addunlock(sched, c->super->grav_up, gravity);
+}
+
+/**
+ * @brief Creates the dependency network for the external gravity tasks of a
+ * given cell.
+ *
+ * @param sched The #scheduler.
+ * @param gravity The gravity task to link.
+ * @param c The cell.
+ */
+static inline void engine_make_external_gravity_dependencies(
+ struct scheduler *sched, struct task *gravity, struct cell *c) {
+
+ /* init --> external gravity --> kick */
+ scheduler_addunlock(sched, c->super->init, gravity);
+ scheduler_addunlock(sched, gravity, c->super->kick);
}
/**
@@ -1474,19 +1455,27 @@ void engine_link_gravity_tasks(struct engine *e) {
/* Gather the multipoles --> mm interaction --> kick */
scheduler_addunlock(sched, gather, t);
- scheduler_addunlock(sched, t, t->ci->gsuper->kick);
+ scheduler_addunlock(sched, t, t->ci->super->kick);
/* init --> mm interaction */
- scheduler_addunlock(sched, t->ci->gsuper->init, t);
+ scheduler_addunlock(sched, t->ci->super->init, t);
}
- /* Self-interaction? */
+ /* Self-interaction for self-gravity? */
if (t->type == task_type_self && t->subtype == task_subtype_grav) {
engine_make_gravity_dependencies(sched, t, t->ci);
}
+ /* Self-interaction for external gravity ? */
+ else if (t->type == task_type_self &&
+ t->subtype == task_subtype_external_grav) {
+
+ engine_make_external_gravity_dependencies(sched, t, t->ci);
+
+ }
+
/* Otherwise, pair interaction? */
else if (t->type == task_type_pair && t->subtype == task_subtype_grav) {
@@ -1495,7 +1484,7 @@ void engine_link_gravity_tasks(struct engine *e) {
engine_make_gravity_dependencies(sched, t, t->ci);
}
- if (t->cj->nodeID == nodeID && t->ci->gsuper != t->cj->gsuper) {
+ if (t->cj->nodeID == nodeID && t->ci->super != t->cj->super) {
engine_make_gravity_dependencies(sched, t, t->cj);
}
@@ -1510,6 +1499,15 @@ void engine_link_gravity_tasks(struct engine *e) {
}
}
+ /* Sub-self-interaction for external gravity ? */
+ else if (t->type == task_type_sub_self &&
+ t->subtype == task_subtype_external_grav) {
+
+ if (t->ci->nodeID == nodeID) {
+ engine_make_external_gravity_dependencies(sched, t, t->ci);
+ }
+ }
+
/* Otherwise, sub-pair interaction? */
else if (t->type == task_type_sub_pair && t->subtype == task_subtype_grav) {
@@ -1517,7 +1515,7 @@ void engine_link_gravity_tasks(struct engine *e) {
engine_make_gravity_dependencies(sched, t, t->ci);
}
- if (t->cj->nodeID == nodeID && t->ci->gsuper != t->cj->gsuper) {
+ if (t->cj->nodeID == nodeID && t->ci->super != t->cj->super) {
engine_make_gravity_dependencies(sched, t, t->cj);
}
@@ -1781,12 +1779,6 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
}
#endif
}
-
- /* External gravity tasks should depend on init and unlock the kick */
- else if (t->type == task_type_grav_external) {
- scheduler_addunlock(sched, t->ci->init, t);
- scheduler_addunlock(sched, t, t->ci->kick);
- }
/* Cooling tasks should depend on kick and unlock sourceterms */
else if (t->type == task_type_cooling) {
scheduler_addunlock(sched, t->ci->kick, t);
@@ -1862,6 +1854,13 @@ void engine_maketasks(struct engine *e) {
/* Add the gravity mm tasks. */
if (e->policy & engine_policy_self_gravity) engine_make_gravity_tasks(e);
+ /* Add the external gravity tasks. */
+ if (e->policy & engine_policy_external_gravity)
+ engine_make_external_gravity_tasks(e);
+
+ if (e->sched.nr_tasks == 0 && (s->nr_gparts > 0 || s->nr_parts > 0))
+ error("We have particles but no hydro or gravity tasks were created.");
+
/* Split the tasks. */
scheduler_splittasks(sched);
@@ -1885,25 +1884,24 @@ void engine_maketasks(struct engine *e) {
/* Count the number of tasks associated with each cell and
store the density tasks in each cell, and make each sort
depend on the sorts of its progeny. */
- if (e->policy & engine_policy_hydro) engine_count_and_link_tasks(e);
+ engine_count_and_link_tasks(e);
- /* Append hierarchical tasks to each cells */
- if (e->policy & engine_policy_hydro)
- for (int k = 0; k < nr_cells; k++)
- engine_make_hydro_hierarchical_tasks(e, &cells[k], NULL);
+ /* Now that the self/pair tasks are at the right level, set the super
+ * pointers. */
+ for (int k = 0; k < nr_cells; k++) cell_set_super(&cells[k], NULL);
- if ((e->policy & engine_policy_self_gravity) ||
- (e->policy & engine_policy_external_gravity))
- for (int k = 0; k < nr_cells; k++)
- engine_make_gravity_hierarchical_tasks(e, &cells[k], NULL);
+ /* Append hierarchical tasks to each cells */
+ for (int k = 0; k < nr_cells; k++)
+ engine_make_hierarchical_tasks(e, &cells[k]);
/* Run through the tasks and make force tasks for each density task.
Each force task depends on the cell ghosts and unlocks the kick task
of its super-cell. */
if (e->policy & engine_policy_hydro) engine_make_extra_hydroloop_tasks(e);
- /* Add the dependencies for the self-gravity stuff */
- if (e->policy & engine_policy_self_gravity) engine_link_gravity_tasks(e);
+ /* Add the dependencies for the gravity stuff */
+ if (e->policy & (engine_policy_self_gravity | engine_policy_external_gravity))
+ engine_link_gravity_tasks(e);
#ifdef WITH_MPI
@@ -2625,7 +2623,10 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* Add the tasks corresponding to external gravity to the masks */
if (e->policy & engine_policy_external_gravity) {
- mask |= 1 << task_type_grav_external;
+ mask |= 1 << task_type_self;
+ mask |= 1 << task_type_sub_self;
+
+ submask |= 1 << task_subtype_external_grav;
}
/* Add MPI tasks if need be */
@@ -2794,7 +2795,11 @@ void engine_step(struct engine *e) {
/* Add the tasks corresponding to external gravity to the masks */
if (e->policy & engine_policy_external_gravity) {
- mask |= 1 << task_type_grav_external;
+
+ mask |= 1 << task_type_self;
+ mask |= 1 << task_type_sub_self;
+
+ submask |= 1 << task_subtype_external_grav;
}
/* Add the tasks corresponding to cooling to the masks */
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index a383a2fdce9e64a5673d0136184368220dc08458..bd970795bdf070a7bd7915cc4f493218dbf319d1 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -524,6 +524,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
/* conserved.energy is NOT the specific energy (u), but the total thermal
energy (u*m) */
p->conserved.energy = u * p->conserved.mass;
+ p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u;
}
/**
@@ -540,4 +541,5 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
p->conserved.energy = gas_internal_energy_from_entropy(p->primitives.rho, S) *
p->conserved.mass;
+ p->primitives.P = gas_pressure_from_entropy(p->primitives.rho, S);
}
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 2c037a05fcba1d9f7a491e97895b57ae0e35eac2..0c69728a9ce6317a8d7ddab945e7aab6e94ee247 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -338,11 +338,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute "grad h" term (note we use rho here and not rho_bar !)*/
const float rho_inv = 1.f / p->rho;
- const float entropy_minus_one_over_gamma = 1.f / p->entropy_one_over_gamma;
const float rho_dh =
1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
- const float pressure_dh = p->density.pressure_dh * rho_inv * p->h *
- hydro_dimension_inv * entropy_minus_one_over_gamma;
+ const float pressure_dh =
+ p->density.pressure_dh * rho_inv * p->h * hydro_dimension_inv;
const float grad_h_term = rho_dh * pressure_dh;
@@ -442,7 +441,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
p->force.h_dt *= p->h * hydro_dimension_inv;
- p->entropy_dt *=
+ p->entropy_dt =
0.5f * gas_entropy_from_internal_energy(p->rho_bar, p->entropy_dt);
}
diff --git a/src/hydro/PressureEntropy/hydro_iact.h b/src/hydro/PressureEntropy/hydro_iact.h
index 18e22233f6f53e488119a58a988ba4c9faf3db36..ce1c38ca69954252dc804af9181b9060a14afcb9 100644
--- a/src/hydro/PressureEntropy/hydro_iact.h
+++ b/src/hydro/PressureEntropy/hydro_iact.h
@@ -255,8 +255,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
- const float sph_term = (S_gamma_j / S_gamma_i - f_i) * P_over_rho2_i * wi_dr +
- (S_gamma_i / S_gamma_j - f_j) * P_over_rho2_j * wj_dr;
+ const float sph_term =
+ (S_gamma_j / S_gamma_i - f_i / S_gamma_i) * P_over_rho2_i * wi_dr +
+ (S_gamma_i / S_gamma_j - f_j / S_gamma_j) * P_over_rho2_j * wj_dr;
/* Eventually got the acceleration */
const float acc = (visc_term + sph_term) * r_inv;
@@ -365,8 +366,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
- const float sph_term = (S_gamma_j / S_gamma_i - f_i) * P_over_rho2_i * wi_dr +
- (S_gamma_i / S_gamma_j - f_j) * P_over_rho2_j * wj_dr;
+ const float sph_term =
+ (S_gamma_j / S_gamma_i - f_i / S_gamma_i) * P_over_rho2_i * wi_dr +
+ (S_gamma_i / S_gamma_j - f_j / S_gamma_j) * P_over_rho2_j * wj_dr;
/* Eventually got the acceleration */
const float acc = (visc_term + sph_term) * r_inv;
diff --git a/src/minmax.h b/src/minmax.h
index 9d92cd71d849dba615fdb05bc342014e0593d989..a53093663c79cf4280d136747663552e49c7f1b2 100644
--- a/src/minmax.h
+++ b/src/minmax.h
@@ -43,4 +43,18 @@
_a > _b ? _a : _b; \
})
+/**
+ * @brief Limits the value of x to be between a and b
+ *
+ * Only wraps once. If x > 2b, the returned value will be larger than b.
+ * Similarly for x < -b.
+ */
+#define box_wrap(x, a, b) \
+ ({ \
+ const __typeof__(x) _x = (x); \
+ const __typeof__(a) _a = (a); \
+ const __typeof__(b) _b = (b); \
+ _x < _a ? (_x + _b) : ((_x > _b) ? (_x - _b) : _x); \
+ })
+
#endif /* SWIFT_MINMAX_H */
diff --git a/src/riemann.h b/src/riemann.h
index d0ae57a640e13c2098708735d6c34de70ebea5b0..69e1a038fc37cb07e1094c8e6736f79b6afa4f35 100644
--- a/src/riemann.h
+++ b/src/riemann.h
@@ -27,18 +27,17 @@
#include "stdio.h"
#include "stdlib.h"
-/* Check that we use an ideal equation of state, since other equations of state
- are not compatible with these Riemann solvers. */
-#ifndef EOS_IDEAL_GAS
-#error Currently there are no Riemann solvers that can handle the requested \
- equation of state. Select an ideal gas equation of state if you want to \
- use this hydro scheme!
-#endif
-
#if defined(RIEMANN_SOLVER_EXACT)
#define RIEMANN_SOLVER_IMPLEMENTATION "Exact Riemann solver (Toro 2009)"
+#if defined(EOS_IDEAL_GAS)
#include "riemann/riemann_exact.h"
+#elif defined(EOS_ISOTHERMAL_GAS)
+#include "riemann/riemann_exact_isothermal.h"
+#else
+#error \
+ "The Exact Riemann solver is incompatible with the selected equation of state!"
+#endif
#elif defined(RIEMANN_SOLVER_TRRS)
diff --git a/src/riemann/riemann_exact_isothermal.h b/src/riemann/riemann_exact_isothermal.h
new file mode 100644
index 0000000000000000000000000000000000000000..8e8a6a4dcce456b55f85148bed7342ad4e651c1b
--- /dev/null
+++ b/src/riemann/riemann_exact_isothermal.h
@@ -0,0 +1,450 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *
+ * 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_RIEMANN_EXACT_ISOTHERMAL_H
+#define SWIFT_RIEMANN_EXACT_ISOTHERMAL_H
+
+#include <float.h>
+#include "adiabatic_index.h"
+#include "minmax.h"
+#include "riemann_vacuum.h"
+
+#define const_isothermal_soundspeed \
+ sqrtf(hydro_gamma_minus_one* const_isothermal_internal_energy)
+
+/**
+ * @brief Relative difference between the middle state velocity and the left or
+ * right state velocity used in the middle state density iteration.
+ *
+ * @param rho Current estimate of the middle state density.
+ * @param W Left or right state vector.
+ * @return Density dependent part of the middle state velocity.
+ */
+__attribute__((always_inline)) INLINE static float riemann_fb(float rho,
+ float* W) {
+ if (rho < W[0]) {
+ return const_isothermal_soundspeed * logf(rho / W[0]);
+ } else {
+ return const_isothermal_soundspeed *
+ (sqrtf(rho / W[0]) - sqrtf(W[0] / rho));
+ }
+}
+
+/**
+ * @brief Derivative w.r.t. rho of the function riemann_fb.
+ *
+ * @param rho Current estimate of the middle state density.
+ * @param W Left or right state vector.
+ * @return Derivative of riemann_fb.
+ */
+__attribute__((always_inline)) INLINE static float riemann_fprimeb(float rho,
+ float* W) {
+ if (rho < W[0]) {
+ return const_isothermal_soundspeed * W[0] / rho;
+ } else {
+ return 0.5 * const_isothermal_soundspeed *
+ (sqrtf(rho / W[0]) + sqrtf(W[0] / rho)) / rho;
+ }
+}
+
+/**
+ * @brief Difference between the left and right middle state velocity estimates.
+ *
+ * Since the middle state velocity takes on a constant value, we want to get
+ * this difference as close to zero as possible.
+ *
+ * @param rho Current estimate of the middle state density.
+ * @param WL Left state vector.
+ * @param WR Right state vector.
+ * @param vL Left state velocity along the interface normal.
+ * @param vR Right state velocity along the interface normal.
+ * @return Difference between the left and right middle state velocity
+ * estimates.
+ */
+__attribute__((always_inline)) INLINE static float riemann_f(
+ float rho, float* WL, float* WR, float vL, float vR) {
+ return riemann_fb(rho, WR) + riemann_fb(rho, WL) + vR - vL;
+}
+
+/**
+ * @brief Derivative of riemann_f w.r.t. rho.
+ *
+ * @param rho Current estimate of the middle state density.
+ * @param WL Left state vector.
+ * @param WR Right state vector.
+ * @return Derivative of riemann_f.
+ */
+__attribute__((always_inline)) INLINE static float riemann_fprime(float rho,
+ float* WL,
+ float* WR) {
+ return riemann_fprimeb(rho, WL) + riemann_fprimeb(rho, WR);
+}
+
+/**
+ * @brief Get a good first guess for the middle state density.
+ *
+ * @param WL The left state vector
+ * @param WR The right state vector
+ * @param vL The left velocity along the interface normal
+ * @param vR The right velocity along the interface normal
+ */
+__attribute__((always_inline)) INLINE static float riemann_guess_rho(float* WL,
+ float* WR,
+ float vL,
+ float vR) {
+
+ /* Currently three possibilities and not really an algorithm to decide which
+ one to choose: */
+ /* just the average */
+ return 0.5f * (WL[0] + WR[0]);
+
+ /* two rarefaction approximation */
+ return sqrtf(WL[0] * WR[0] * expf((vL - vR) / const_isothermal_soundspeed));
+
+ /* linearized primitive variable approximation */
+ return 0.25f * (WL[0] + WR[0]) * (vL - vR) / const_isothermal_soundspeed +
+ 0.5f * (WL[0] + WR[0]);
+}
+
+/**
+ * @brief Find the zeropoint of riemann_f(rho) using Brent's method.
+ *
+ * @param lower_limit Lower limit for the method (riemann_f(lower_limit) < 0)
+ * @param upper_limit Upper limit for the method (riemann_f(upper_limit) > 0)
+ * @param lowf Value of riemann_f(lower_limit).
+ * @param upf Value of riemann_f(upper_limit).
+ * @param error_tol Tolerance used to decide if the solution is converged
+ * @param WL Left state vector
+ * @param WR Right state vector
+ * @param vL The left velocity along the interface normal
+ * @param vR The right velocity along the interface normal
+ */
+__attribute__((always_inline)) INLINE static float riemann_solve_brent(
+ float lower_limit, float upper_limit, float lowf, float upf,
+ float error_tol, float* WL, float* WR, float vL, float vR) {
+
+ float a, b, c, d, s;
+ float fa, fb, fc, fs;
+ float tmp, tmp2;
+ int mflag;
+ int i;
+
+ a = lower_limit;
+ b = upper_limit;
+ c = 0.0f;
+ d = FLT_MAX;
+
+ fa = lowf;
+ fb = upf;
+
+ fc = 0.0f;
+ s = 0.0f;
+ fs = 0.0f;
+
+ /* if f(a) f(b) >= 0 then error-exit */
+ if (fa * fb >= 0.0f) {
+ error(
+ "Brent's method called with equal sign function values!\n"
+ "f(%g) = %g, f(%g) = %g\n",
+ a, fa, b, fb);
+ /* return NaN */
+ return 0.0f / 0.0f;
+ }
+
+ /* if |f(a)| < |f(b)| then swap (a,b) */
+ if (fabs(fa) < fabs(fb)) {
+ tmp = a;
+ a = b;
+ b = tmp;
+ tmp = fa;
+ fa = fb;
+ fb = tmp;
+ }
+
+ c = a;
+ fc = fa;
+ mflag = 1;
+ i = 0;
+
+ while (!(fb == 0.0f) && (fabs(a - b) > error_tol * 0.5f * (a + b))) {
+ if ((fa != fc) && (fb != fc)) /* Inverse quadratic interpolation */
+ s = a * fb * fc / (fa - fb) / (fa - fc) +
+ b * fa * fc / (fb - fa) / (fb - fc) +
+ c * fa * fb / (fc - fa) / (fc - fb);
+ else
+ /* Secant Rule */
+ s = b - fb * (b - a) / (fb - fa);
+
+ tmp2 = 0.25f * (3.0f * a + b);
+ if (!(((s > tmp2) && (s < b)) || ((s < tmp2) && (s > b))) ||
+ (mflag && (fabs(s - b) >= (0.5f * fabs(b - c)))) ||
+ (!mflag && (fabs(s - b) >= (0.5f * fabs(c - d)))) ||
+ (mflag && (fabs(b - c) < error_tol)) ||
+ (!mflag && (fabs(c - d) < error_tol))) {
+ s = 0.5f * (a + b);
+ mflag = 1;
+ } else {
+ mflag = 0;
+ }
+ fs = riemann_f(s, WL, WR, vL, vR);
+ d = c;
+ c = b;
+ fc = fb;
+ if (fa * fs < 0.) {
+ b = s;
+ fb = fs;
+ } else {
+ a = s;
+ fa = fs;
+ }
+
+ /* if |f(a)| < |f(b)| then swap (a,b) */
+ if (fabs(fa) < fabs(fb)) {
+ tmp = a;
+ a = b;
+ b = tmp;
+ tmp = fa;
+ fa = fb;
+ fb = tmp;
+ }
+ i++;
+ }
+ return b;
+}
+
+/**
+ * @brief Solve the Riemann problem between the given left and right state and
+ * along the given interface normal
+ *
+ * @param WL The left state vector
+ * @param WR The right state vector
+ * @param Whalf Empty state vector in which the result will be stored
+ * @param n_unit Normal vector of the interface
+ */
+__attribute__((always_inline)) INLINE static void riemann_solver_solve(
+ float* WL, float* WR, float* Whalf, float* n_unit) {
+
+ /* velocity of the left and right state in a frame aligned with n_unit */
+ float vL, vR, vhalf;
+ /* variables used for finding rhostar */
+ float rho, rhoguess, frho, frhoguess;
+ /* variables used for sampling the solution */
+ float u, S, SH, ST;
+
+ int errorFlag = 0;
+
+ /* sanity checks */
+ if (WL[0] != WL[0]) {
+ printf("NaN WL!\n");
+ errorFlag = 1;
+ }
+ if (WR[0] != WR[0]) {
+ printf("NaN WR!\n");
+ errorFlag = 1;
+ }
+ if (WL[0] < 0.0f) {
+ printf("Negative WL!\n");
+ errorFlag = 1;
+ }
+ if (WR[0] < 0.0f) {
+ printf("Negative WR!\n");
+ errorFlag = 1;
+ }
+ if (errorFlag) {
+ printf("WL: %g %g %g %g %g\n", WL[0], WL[1], WL[2], WL[3], WL[4]);
+ printf("WR: %g %g %g %g %g\n", WR[0], WR[1], WR[2], WR[3], WR[4]);
+ error("Riemman solver input error!\n");
+ }
+
+ /* calculate velocities in interface frame */
+ vL = WL[1] * n_unit[0] + WL[2] * n_unit[1] + WL[3] * n_unit[2];
+ vR = WR[1] * n_unit[0] + WR[2] * n_unit[1] + WR[3] * n_unit[2];
+
+ /* VACUUM... */
+
+ rho = 0.;
+ /* obtain a first guess for p */
+ rhoguess = riemann_guess_rho(WL, WR, vL, vR);
+ frho = riemann_f(rho, WL, WR, vL, vR);
+ frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
+ /* ok, rhostar is close to 0, better use Brent's method... */
+ /* we use Newton-Raphson until we find a suitable interval */
+ if (frho * frhoguess >= 0.0f) {
+ /* Newton-Raphson until convergence or until suitable interval is found
+ to use Brent's method */
+ unsigned int counter = 0;
+ while (fabs(rho - rhoguess) > 1.e-6f * 0.5f * (rho + rhoguess) &&
+ frhoguess < 0.0f) {
+ rho = rhoguess;
+ rhoguess = rhoguess - frhoguess / riemann_fprime(rhoguess, WL, WR);
+ frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
+ counter++;
+ if (counter > 1000) {
+ error("Stuck in Newton-Raphson!\n");
+ }
+ }
+ }
+ /* As soon as there is a suitable interval: use Brent's method */
+ if (1.e6 * fabs(rho - rhoguess) > 0.5f * (rho + rhoguess) &&
+ frhoguess > 0.0f) {
+ rho = 0.0f;
+ frho = riemann_f(rho, WL, WR, vL, vR);
+ /* use Brent's method to find the zeropoint */
+ rho = riemann_solve_brent(rho, rhoguess, frho, frhoguess, 1.e-6, WL, WR, vL,
+ vR);
+ } else {
+ rho = rhoguess;
+ }
+
+ /* calculate the middle state velocity */
+ u = 0.5f * (vL - riemann_fb(rho, WL) + vR + riemann_fb(rho, WR));
+
+ /* sample the solution */
+ if (u > 0.0f) {
+ /* left state */
+ Whalf[1] = WL[1];
+ Whalf[2] = WL[2];
+ Whalf[3] = WL[3];
+ if (WL[0] < rho) {
+ /* left shock wave */
+ S = vL - const_isothermal_soundspeed * sqrtf(rho / WL[0]);
+ if (S >= 0.) {
+ /* to the left of the shock */
+ Whalf[0] = WL[0];
+ vhalf = 0.0f;
+ } else {
+ /* to the right of the shock */
+ Whalf[0] = rho;
+ vhalf = u - vL;
+ }
+ } else {
+ /* left rarefaction wave */
+ SH = vL - const_isothermal_soundspeed;
+ ST = u - const_isothermal_soundspeed;
+ if (SH > 0.) {
+ /* to the left of the rarefaction */
+ Whalf[0] = WL[0];
+ vhalf = 0.0f;
+ } else if (ST > 0.0f) {
+ /* inside the rarefaction */
+ Whalf[0] = WL[0] * expf(vL / const_isothermal_soundspeed - 1.0f);
+ vhalf = const_isothermal_soundspeed - vL;
+ } else {
+ /* to the right of the rarefaction */
+ Whalf[0] = rho;
+ vhalf = u - vL;
+ }
+ }
+ } else {
+ /* right state */
+ Whalf[1] = WR[1];
+ Whalf[2] = WR[2];
+ Whalf[3] = WR[3];
+ if (WR[0] < rho) {
+ /* right shock wave */
+ S = vR + const_isothermal_soundspeed * sqrtf(rho / WR[0]);
+ if (S > 0.0f) {
+ /* to the left of the shock wave: middle state */
+ Whalf[0] = rho;
+ vhalf = u - vR;
+ } else {
+ /* to the right of the shock wave: right state */
+ Whalf[0] = WR[0];
+ vhalf = 0.0f;
+ }
+ } else {
+ /* right rarefaction wave */
+ SH = vR + const_isothermal_soundspeed;
+ ST = u + const_isothermal_soundspeed;
+ if (ST > 0.0f) {
+ /* to the left of rarefaction: middle state */
+ Whalf[0] = rho;
+ vhalf = u - vR;
+ } else if (SH > 0.0f) {
+ /* inside rarefaction */
+ Whalf[0] = WR[0] * expf(-vR / const_isothermal_soundspeed - 1.0f);
+ vhalf = -const_isothermal_soundspeed - vR;
+ } else {
+ /* to the right of rarefaction: right state */
+ Whalf[0] = WR[0];
+ vhalf = 0.0f;
+ }
+ }
+ }
+
+ /* add the velocity solution along the interface normal to the velocities */
+ Whalf[1] += vhalf * n_unit[0];
+ Whalf[2] += vhalf * n_unit[1];
+ Whalf[3] += vhalf * n_unit[2];
+
+ /* the pressure is completely irrelevant in this case */
+ Whalf[4] =
+ Whalf[0] * const_isothermal_soundspeed * const_isothermal_soundspeed;
+}
+
+__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
+ float* Wi, float* Wj, float* n_unit, float* vij, float* totflux) {
+
+ float Whalf[5];
+ float flux[5][3];
+ float vtot[3];
+ float rhoe;
+
+ riemann_solver_solve(Wi, Wj, Whalf, n_unit);
+
+ flux[0][0] = Whalf[0] * Whalf[1];
+ flux[0][1] = Whalf[0] * Whalf[2];
+ flux[0][2] = Whalf[0] * Whalf[3];
+
+ vtot[0] = Whalf[1] + vij[0];
+ vtot[1] = Whalf[2] + vij[1];
+ vtot[2] = Whalf[3] + vij[2];
+ flux[1][0] = Whalf[0] * vtot[0] * Whalf[1] + Whalf[4];
+ flux[1][1] = Whalf[0] * vtot[0] * Whalf[2];
+ flux[1][2] = Whalf[0] * vtot[0] * Whalf[3];
+ flux[2][0] = Whalf[0] * vtot[1] * Whalf[1];
+ flux[2][1] = Whalf[0] * vtot[1] * Whalf[2] + Whalf[4];
+ flux[2][2] = Whalf[0] * vtot[1] * Whalf[3];
+ flux[3][0] = Whalf[0] * vtot[2] * Whalf[1];
+ flux[3][1] = Whalf[0] * vtot[2] * Whalf[2];
+ flux[3][2] = Whalf[0] * vtot[2] * Whalf[3] + Whalf[4];
+
+ /* eqn. (15) */
+ /* F_P = \rho e ( \vec{v} - \vec{v_{ij}} ) + P \vec{v} */
+ /* \rho e = P / (\gamma-1) + 1/2 \rho \vec{v}^2 */
+ rhoe = Whalf[4] / hydro_gamma_minus_one +
+ 0.5f * Whalf[0] *
+ (vtot[0] * vtot[0] + vtot[1] * vtot[1] + vtot[2] * vtot[2]);
+ flux[4][0] = rhoe * Whalf[1] + Whalf[4] * vtot[0];
+ flux[4][1] = rhoe * Whalf[2] + Whalf[4] * vtot[1];
+ flux[4][2] = rhoe * Whalf[3] + Whalf[4] * vtot[2];
+
+ totflux[0] =
+ flux[0][0] * n_unit[0] + flux[0][1] * n_unit[1] + flux[0][2] * n_unit[2];
+ totflux[1] =
+ flux[1][0] * n_unit[0] + flux[1][1] * n_unit[1] + flux[1][2] * n_unit[2];
+ totflux[2] =
+ flux[2][0] * n_unit[0] + flux[2][1] * n_unit[1] + flux[2][2] * n_unit[2];
+ totflux[3] =
+ flux[3][0] * n_unit[0] + flux[3][1] * n_unit[1] + flux[3][2] * n_unit[2];
+ totflux[4] =
+ flux[4][0] * n_unit[0] + flux[4][1] * n_unit[1] + flux[4][2] * n_unit[2];
+}
+
+#endif /* SWIFT_RIEMANN_EXACT_ISOTHERMAL_H */
diff --git a/src/riemann/riemann_hllc.h b/src/riemann/riemann_hllc.h
index 1a4c3f5f8338f6504a8c5f1d9eab8451eb20246c..962fb8380ad9919d7b7ad96325a12f0bfc53d255 100644
--- a/src/riemann/riemann_hllc.h
+++ b/src/riemann/riemann_hllc.h
@@ -24,6 +24,11 @@
#include "minmax.h"
#include "riemann_vacuum.h"
+#ifndef EOS_IDEAL_GAS
+#error \
+ "The HLLC Riemann solver currently only supports and ideal gas equation of state. Either select this equation of state, or try using another Riemann solver!"
+#endif
+
__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
float *WL, float *WR, float *n, float *vij, float *totflux) {
diff --git a/src/riemann/riemann_trrs.h b/src/riemann/riemann_trrs.h
index b13a76b4c57af548497780e974e5c9ee3a721fac..be56b93583b693d7d492e0c8b1357cdbe57e88b5 100644
--- a/src/riemann/riemann_trrs.h
+++ b/src/riemann/riemann_trrs.h
@@ -23,6 +23,11 @@
#include "adiabatic_index.h"
#include "riemann_vacuum.h"
+#ifndef EOS_IDEAL_GAS
+#error \
+ "The TRRS Riemann solver currently only supports an ideal gas equation of state. Either select this equation of state, or try using another Riemann solver!"
+#endif
+
/**
* @brief Solve the Riemann problem using the Two Rarefaction Riemann Solver
*
diff --git a/src/runner.c b/src/runner.c
index 7294d0e1c862af4f31c0ef32b5f936df56e9df70..cf70e482bc5e31d35500ec00a49e4d888d6270f5 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -184,12 +184,12 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
for (int i = 0; i < gcount; i++) {
/* Get a direct pointer on the part. */
- struct gpart *restrict g = &gparts[i];
+ struct gpart *restrict gp = &gparts[i];
/* Is this part within the time step? */
- if (g->ti_end <= ti_current) {
+ if (gp->ti_end <= ti_current) {
- external_gravity_acceleration(time, potential, constants, g);
+ external_gravity_acceleration(time, potential, constants, gp);
}
}
@@ -1333,9 +1333,12 @@ void *runner_main(void *data) {
runner_doself2_force(r, ci);
else if (t->subtype == task_subtype_grav)
runner_doself_grav(r, ci, 1);
+ else if (t->subtype == task_subtype_external_grav)
+ runner_do_grav_external(r, ci, 1);
else
error("Unknown task subtype.");
break;
+
case task_type_pair:
if (t->subtype == task_subtype_density)
runner_dopair1_density(r, ci, cj);
@@ -1350,9 +1353,7 @@ void *runner_main(void *data) {
else
error("Unknown task subtype.");
break;
- case task_type_sort:
- runner_do_sort(r, ci, t->flags, 1);
- break;
+
case task_type_sub_self:
if (t->subtype == task_subtype_density)
runner_dosub_self1_density(r, ci, 1);
@@ -1364,9 +1365,12 @@ void *runner_main(void *data) {
runner_dosub_self2_force(r, ci, 1);
else if (t->subtype == task_subtype_grav)
runner_dosub_grav(r, ci, cj, 1);
+ else if (t->subtype == task_subtype_external_grav)
+ runner_do_grav_external(r, ci, 1);
else
error("Unknown task subtype.");
break;
+
case task_type_sub_pair:
if (t->subtype == task_subtype_density)
runner_dosub_pair1_density(r, ci, cj, t->flags, 1);
@@ -1381,6 +1385,10 @@ void *runner_main(void *data) {
else
error("Unknown task subtype.");
break;
+
+ case task_type_sort:
+ runner_do_sort(r, ci, t->flags, 1);
+ break;
case task_type_init:
runner_do_init(r, ci, 1);
break;
@@ -1424,9 +1432,6 @@ void *runner_main(void *data) {
case task_type_grav_fft:
runner_do_grav_fft(r);
break;
- case task_type_grav_external:
- runner_do_grav_external(r, t->ci, 1);
- break;
case task_type_cooling:
runner_do_cooling(r, t->ci, 1);
break;
diff --git a/src/scheduler.c b/src/scheduler.c
index b4e2575807aff4cb8b3878dcb25a773269038b74..e3fe92e51ad810ad932cdd307dd25c7fe6f60f4c 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1185,20 +1185,10 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
break;
case task_type_pair:
case task_type_sub_pair:
- if (t->subtype == task_subtype_grav) {
-
- qid = t->ci->gsuper->owner;
- if (qid < 0 ||
- s->queues[qid].count > s->queues[t->cj->gsuper->owner].count)
- qid = t->cj->gsuper->owner;
-
- } else {
-
- qid = t->ci->super->owner;
- if (qid < 0 ||
- s->queues[qid].count > s->queues[t->cj->super->owner].count)
- qid = t->cj->super->owner;
- }
+ qid = t->ci->super->owner;
+ if (qid < 0 ||
+ s->queues[qid].count > s->queues[t->cj->super->owner].count)
+ qid = t->cj->super->owner;
break;
case task_type_recv:
#ifdef WITH_MPI
diff --git a/src/space.c b/src/space.c
index ddd4c076e9ef9c18c8bf25eeb9360dc7333b88f8..186e11c066f3c8d788c3eaa189da3ed2461aed87 100644
--- a/src/space.c
+++ b/src/space.c
@@ -112,6 +112,15 @@ struct parallel_sort {
volatile unsigned int first, last, waiting;
};
+/**
+ * @brief Information required to compute the particle cell indices.
+ */
+struct index_data {
+ struct space *s;
+ struct cell *cells;
+ int *ind;
+};
+
/**
* @brief Get the shift-id of the given pair of cells, swapping them
* if need be.
@@ -326,7 +335,6 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
c->count = 0;
c->gcount = 0;
c->super = c;
- c->gsuper = c;
c->ti_old = ti_current;
lock_init(&c->lock);
}
@@ -402,7 +410,6 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
s->cells_top[k].cooling = NULL;
s->cells_top[k].sourceterms = NULL;
s->cells_top[k].super = &s->cells_top[k];
- s->cells_top[k].gsuper = &s->cells_top[k];
#if WITH_MPI
s->cells_top[k].recv_xv = NULL;
s->cells_top[k].recv_rho = NULL;
@@ -446,49 +453,21 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
struct cell *restrict cells_top = s->cells_top;
const int ti_current = (s->e != NULL) ? s->e->ti_current : 0;
- const double ih[3] = {s->iwidth[0], s->iwidth[1], s->iwidth[2]};
- const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
- const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
-
/* Run through the particles and get their cell index. */
- // tic = getticks();
const size_t ind_size = s->size_parts;
int *ind;
if ((ind = (int *)malloc(sizeof(int) * ind_size)) == NULL)
error("Failed to allocate temporary particle indices.");
- for (size_t k = 0; k < nr_parts; k++) {
- struct part *restrict p = &s->parts[k];
- for (int j = 0; j < 3; j++)
- if (p->x[j] < 0.0)
- p->x[j] += dim[j];
- else if (p->x[j] >= dim[j])
- p->x[j] -= dim[j];
- ind[k] =
- cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
- cells_top[ind[k]].count++;
- }
- // message( "getting particle indices took %.3f %s." ,
- // clocks_from_ticks(getticks() - tic), clocks_getunit()):
+ if (ind_size > 0) space_parts_get_cell_index(s, ind, cells_top, verbose);
+ for (size_t i = 0; i < s->nr_parts; ++i) cells_top[ind[i]].count++;
/* Run through the gravity particles and get their cell index. */
- // tic = getticks();
const size_t gind_size = s->size_gparts;
int *gind;
if ((gind = (int *)malloc(sizeof(int) * gind_size)) == NULL)
error("Failed to allocate temporary g-particle indices.");
- for (size_t k = 0; k < nr_gparts; k++) {
- struct gpart *restrict gp = &s->gparts[k];
- for (int j = 0; j < 3; j++)
- if (gp->x[j] < 0.0)
- gp->x[j] += dim[j];
- else if (gp->x[j] >= dim[j])
- gp->x[j] -= dim[j];
- gind[k] =
- cell_getid(cdim, gp->x[0] * ih[0], gp->x[1] * ih[1], gp->x[2] * ih[2]);
- cells_top[gind[k]].gcount++;
- }
-// message( "getting g-particle indices took %.3f %s." ,
-// clocks_from_ticks(getticks() - tic), clocks_getunit());
+ if (gind_size > 0) space_gparts_get_cell_index(s, gind, cells_top, verbose);
+ for (size_t i = 0; i < s->nr_gparts; ++i) cells_top[gind[i]].gcount++;
#ifdef WITH_MPI
@@ -592,6 +571,9 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
ind = ind_new;
}
+ const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
+ const double ih[3] = {s->iwidth[0], s->iwidth[1], s->iwidth[2]};
+
/* Assign each particle to its cell. */
for (size_t k = nr_parts; k < s->nr_parts; k++) {
const struct part *const p = &s->parts[k];
@@ -619,9 +601,10 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
for (size_t k = 1; k < nr_parts; k++) {
if (ind[k - 1] > ind[k]) {
error("Sort failed!");
- } else if (ind[k] != cell_getid(cdim, s->parts[k].x[0] * ih[0],
- s->parts[k].x[1] * ih[1],
- s->parts[k].x[2] * ih[2])) {
+ } else if (ind[k] != cell_getid(s->cdim,
+ s->parts[k].x[0] * s->iwidth[0],
+ s->parts[k].x[1] * s->iwidth[1],
+ s->parts[k].x[2] * s->iwidth[2])) {
error("Incorrect indices!");
}
}
@@ -766,6 +749,164 @@ void space_sanitize(struct space *s) {
}
}
+/**
+ * @brief #threadpool mapper function to compute the particle cell indices.
+ *
+ * @param map_data Pointer towards the particles.
+ * @param nr_parts The number of particles to treat.
+ * @param extra_data Pointers to the space and index list
+ */
+void space_parts_get_cell_index_mapper(void *map_data, int nr_parts,
+ void *extra_data) {
+
+ /* Unpack the data */
+ struct part *restrict parts = (struct part *)map_data;
+ struct index_data *data = (struct index_data *)extra_data;
+ struct space *s = data->s;
+ int *const ind = data->ind + (ptrdiff_t)(parts - s->parts);
+
+ /* Get some constants */
+ const double dim_x = s->dim[0];
+ const double dim_y = s->dim[1];
+ const double dim_z = s->dim[2];
+ const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
+ const double ih_x = s->iwidth[0];
+ const double ih_y = s->iwidth[1];
+ const double ih_z = s->iwidth[2];
+
+ for (int k = 0; k < nr_parts; k++) {
+
+ /* Get the particle */
+ struct part *restrict p = &parts[k];
+
+ const double old_pos_x = p->x[0];
+ const double old_pos_y = p->x[1];
+ const double old_pos_z = p->x[2];
+
+ /* Put it back into the simulation volume */
+ const double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
+ const double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
+ const double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+
+ /* Get its cell index */
+ const int index =
+ cell_getid(cdim, pos_x * ih_x, pos_y * ih_y, pos_z * ih_z);
+ ind[k] = index;
+
+ /* Update the position */
+ p->x[0] = pos_x;
+ p->x[1] = pos_y;
+ p->x[2] = pos_z;
+ }
+}
+
+/**
+ * @brief #threadpool mapper function to compute the g-particle cell indices.
+ *
+ * @param map_data Pointer towards the g-particles.
+ * @param nr_gparts The number of g-particles to treat.
+ * @param extra_data Pointers to the space and index list
+ */
+void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts,
+ void *extra_data) {
+
+ /* Unpack the data */
+ struct gpart *restrict gparts = (struct gpart *)map_data;
+ struct index_data *data = (struct index_data *)extra_data;
+ struct space *s = data->s;
+ int *const ind = data->ind + (ptrdiff_t)(gparts - s->gparts);
+
+ /* Get some constants */
+ const double dim_x = s->dim[0];
+ const double dim_y = s->dim[1];
+ const double dim_z = s->dim[2];
+ const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
+ const double ih_x = s->iwidth[0];
+ const double ih_y = s->iwidth[1];
+ const double ih_z = s->iwidth[2];
+
+ for (int k = 0; k < nr_gparts; k++) {
+
+ /* Get the particle */
+ struct gpart *restrict gp = &gparts[k];
+
+ const double old_pos_x = gp->x[0];
+ const double old_pos_y = gp->x[1];
+ const double old_pos_z = gp->x[2];
+
+ /* Put it back into the simulation volume */
+ const double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
+ const double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
+ const double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+
+ /* Get its cell index */
+ const int index =
+ cell_getid(cdim, pos_x * ih_x, pos_y * ih_y, pos_z * ih_z);
+ ind[k] = index;
+
+ /* Update the position */
+ gp->x[0] = pos_x;
+ gp->x[1] = pos_y;
+ gp->x[2] = pos_z;
+ }
+}
+
+/**
+ * @brief Computes the cell index of all the particles and update the cell
+ * count.
+ *
+ * @param s The #space.
+ * @param ind The array of indices to fill.
+ * @param cells The array of #cell to update.
+ * @param verbose Are we talkative ?
+ */
+void space_parts_get_cell_index(struct space *s, int *ind, struct cell *cells,
+ int verbose) {
+
+ const ticks tic = getticks();
+
+ /* Pack the extra information */
+ struct index_data data;
+ data.s = s;
+ data.cells = cells;
+ data.ind = ind;
+
+ threadpool_map(&s->e->threadpool, space_parts_get_cell_index_mapper, s->parts,
+ s->nr_parts, sizeof(struct part), 1000, &data);
+
+ if (verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
+
+/**
+ * @brief Computes the cell index of all the g-particles and update the cell
+ * gcount.
+ *
+ * @param s The #space.
+ * @param gind The array of indices to fill.
+ * @param cells The array of #cell to update.
+ * @param verbose Are we talkative ?
+ */
+void space_gparts_get_cell_index(struct space *s, int *gind, struct cell *cells,
+ int verbose) {
+
+ const ticks tic = getticks();
+
+ /* Pack the extra information */
+ struct index_data data;
+ data.s = s;
+ data.cells = cells;
+ data.ind = gind;
+
+ threadpool_map(&s->e->threadpool, space_gparts_get_cell_index_mapper,
+ s->gparts, s->nr_gparts, sizeof(struct gpart), 1000, &data);
+
+ if (verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
+
/**
* @brief Sort the particles and condensed particles according to the given
* indices.
@@ -1351,7 +1492,6 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
temp->nodeID = c->nodeID;
temp->parent = c;
temp->super = NULL;
- temp->gsuper = NULL;
c->progeny[k] = temp;
}
diff --git a/src/space.h b/src/space.h
index 85f32bd7580945de6d9713316b557c92667987c0..011dfb71a6c3ac2b51093ce83bc6b65ceecc2821 100644
--- a/src/space.h
+++ b/src/space.h
@@ -170,6 +170,10 @@ void space_recycle(struct space *s, struct cell *c);
void space_split(struct space *s, struct cell *cells, int nr_cells,
int verbose);
void space_split_mapper(void *map_data, int num_elements, void *extra_data);
+void space_parts_get_cell_index(struct space *s, int *ind, struct cell *cells,
+ int verbose);
+void space_gparts_get_cell_index(struct space *s, int *gind, struct cell *cells,
+ int verbose);
void space_do_parts_sort();
void space_do_gparts_sort();
void space_init_parts(struct space *s);
diff --git a/src/task.c b/src/task.c
index 00068f45769b4ca606cc729bd5e89c13ae729eef..54b9363b7ac3d5c372b591c00b0b03cc274f66b5 100644
--- a/src/task.c
+++ b/src/task.c
@@ -48,13 +48,13 @@
/* Task type names. */
const char *taskID_names[task_type_count] = {
- "none", "sort", "self", "pair", "sub_self",
- "sub_pair", "init", "ghost", "extra_ghost", "kick",
- "kick_fixdt", "send", "recv", "grav_gather_m", "grav_fft",
- "grav_mm", "grav_up", "grav_external", "cooling", "sourceterms"};
+ "none", "sort", "self", "pair", "sub_self",
+ "sub_pair", "init", "ghost", "extra_ghost", "kick",
+ "kick_fixdt", "send", "recv", "grav_gather_m", "grav_fft",
+ "grav_mm", "grav_up", "cooling", "sourceterms"};
const char *subtaskID_names[task_subtype_count] = {
- "none", "density", "gradient", "force", "grav", "tend"};
+ "none", "density", "gradient", "force", "grav", "external_grav", "tend"};
/**
* @brief Computes the overlap between the parts array of two given cells.
@@ -135,6 +135,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
break;
case task_subtype_grav:
+ case task_subtype_external_grav:
return task_action_gpart;
break;
@@ -150,7 +151,14 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_type_kick_fixdt:
case task_type_send:
case task_type_recv:
- return task_action_all;
+ if (t->ci->count > 0 && t->ci->gcount > 0)
+ return task_action_all;
+ else if (t->ci->count > 0)
+ return task_action_part;
+ else if (t->ci->gcount > 0)
+ return task_action_gpart;
+ else
+ error("Task without particles");
break;
case task_type_grav_gather_m:
@@ -160,10 +168,6 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
return task_action_multipole;
break;
- case task_type_grav_external:
- return task_action_gpart;
- break;
-
default:
error("Unknown task_action for task");
return task_action_none;
@@ -393,3 +397,15 @@ void task_print_submask(unsigned int submask) {
printf(" %s=%s", subtaskID_names[k], (submask & (1 << k)) ? "yes" : "no");
printf(" ]\n");
}
+
+/**
+ * @brief Print basic information about a task.
+ *
+ * @param t The #task.
+ */
+void task_print(const struct task *t) {
+
+ message("Type:'%s' sub_type:'%s' wait=%d nr_unlocks=%d skip=%d",
+ taskID_names[t->type], subtaskID_names[t->subtype], t->wait,
+ t->nr_unlock_tasks, t->skip);
+}
diff --git a/src/task.h b/src/task.h
index bc4df3dc2a4cfee3c382e9f2059cba84f29299f7..f840c0b4b8e807dce28f6f13479dbdf4995ab66d 100644
--- a/src/task.h
+++ b/src/task.h
@@ -53,7 +53,6 @@ enum task_types {
task_type_grav_fft,
task_type_grav_mm,
task_type_grav_up,
- task_type_grav_external,
task_type_cooling,
task_type_sourceterms,
task_type_count
@@ -68,6 +67,7 @@ enum task_subtypes {
task_subtype_gradient,
task_subtype_force,
task_subtype_grav,
+ task_subtype_external_grav,
task_subtype_tend,
task_subtype_count
} __attribute__((packed));
@@ -160,5 +160,6 @@ int task_lock(struct task *t);
void task_print_mask(unsigned int mask);
void task_print_submask(unsigned int submask);
void task_do_rewait(struct task *t);
+void task_print(const struct task *t);
#endif /* SWIFT_TASK_H */
diff --git a/theory/SPH/Flavours/sph_flavours.tex b/theory/SPH/Flavours/sph_flavours.tex
index 1ca7fb7018ed14cbf4972a181b66f841b75b6f43..5fe1277373552d60607671299437a371e068169c 100644
--- a/theory/SPH/Flavours/sph_flavours.tex
+++ b/theory/SPH/Flavours/sph_flavours.tex
@@ -305,17 +305,72 @@ P_{\partial h_i}$ and $\rho_{\partial h_i}$
(eq. \ref{eq:sph:minimal:rho_dh}):
\begin{equation}
- f_i \equiv \left(\frac{\bar P_{\partial h_i} h_i}{3\rho_i
- \tilde{A}_i}\right)\left(1 + \frac{h_i}{3\rho_i}\rho_{\partial
- h_i}\right)^{-1}.
+ f_i \equiv \left(\frac{h_i}{3\rho_i}\bar P_{\partial
+ h_i}\right)\left(1 + \frac{h_i}{3\rho_i}\rho_{\partial
+ h_i}\right)^{-1}.
\end{equation}
\subsubsection{Hydrodynamical accelerations (\nth{2} neighbour loop)}
+The accelerations are given by the following term:
+
+\begin{align}
+ \frac{d\vec{v}_i}{dt} = -\sum_j m_j &\left[\frac{\bar P_i}{\bar\rho_i^2} \left(\frac{\tilde A_j}{\tilde A_i} - \frac{f_i}{\tilde A_i}\right)\nabla_x W(\vec{x}_{ij}, h_i) \right. \nonumber \\
+ &+\frac{P_j}{\rho_j^2} \left(\frac{\tilde A_i}{\tilde A_j} - \frac{f_j}{\tilde A_j}\right)\nabla_x W(\vec{x}_{ij},h_j) \\
+ &+ \left. \bigg.\nu_{ij} \Wij \right], \label{eq:sph:pe:dv_dt}
+\end{align}
+where the viscosity term $\nu_{ij}$ has been computed as in
+the \GadgetSPH case (Eq. \ref{eq:sph:gadget2:balsara}
+and \ref{eq:sph:gadget2:nu_ij}). For completeness, the equation of
+motion for the entropy is
+
+\begin{equation}
+\frac{dA_i}{dt} = \frac{1}{2} A_{\rm eos}\left(\rho_i, \sum_j
+m_j \nu_{ij}\vec{v}_{ij}\cdot \Wij\right).
+\end{equation}
+
+\subsubsection{Time integration}
+
+The time-step condition is identical to the \MinimalSPH case
+(Eq. \ref{eq:sph:minimal:dt}). The same applies to the integration
+forward in time (Eq. \ref{eq:sph:minimal:kick_v} to
+\ref{eq:sph:minimal:kick_c}) with the exception of the change in
+internal energy (Eq. \ref{eq:sph:minimal:kick_u}) which gets replaced
+by an integration for the the entropy:
+
+\begin{align}
+ \vec{v}_i &\rightarrow \vec{v}_i + \frac{d\vec{v}_i}{dt} \Delta t \label{eq:sph:pe:kick_v}\\
+ A_i &\rightarrow A_i + \frac{dA_i}{dt} \Delta t \label{eq:sph:pe:kick_A}\\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i\right) \label{eq:sph:pe:kick_P}, \\
+ c_i &\rightarrow c_{\rm eos}\left(\rho_i,
+ A_i\right) \label{eq:sph:pe:kick_c}, \\
+ \tilde A_i &= A_i^{1/\gamma}
+\end{align}
+where, once again, we made use of the equation of state relating
+thermodynamical quantities.
\subsubsection{Particle properties prediction}
+The prediction step is also identical to the \MinimalSPH case with the
+entropic function replacing the thermal energy.
+
+\begin{align}
+ \vec{x}_i &\rightarrow \vec{x}_i + \vec{v}_i \Delta t \label{eq:sph:pe:drift_x} \\
+ h_i &\rightarrow h_i \exp\left(\frac{1}{h_i} \frac{dh_i}{dt}
+ \Delta t\right), \label{eq:sph:pe:drift_h}\\
+ \rho_i &\rightarrow \rho_i \exp\left(-\frac{3}{h_i} \frac{dh_i}{dt}
+ \Delta t\right), \label{eq:sph:pe:drift_rho} \\
+ \tilde A_i &\rightarrow \left(A_i + \frac{dA_i}{dt}
+ \Delta t\right)^{1/\gamma} \label{eq:sph:pe:drift_A_tilde}, \\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt} \Delta t\right), \label{eq:sph:pe:drift_P}\\
+ c_i &\rightarrow c_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt}
+ \Delta t\right) \label{eq:sph:pe:drift_c},
+\end{align}
+where, as above, the last two updated quantities are obtained using
+the pre-defined equation of state. Note that the entropic function $A_i$
+itself is \emph{not} updated.
+
\subsection{Pressure-Energy SPH}
\label{sec:sph:pu}