Skip to content
GitLab
Commit 1c1d49f4 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge remote-tracking branch 'origin/master' into mark_tasks_in_drift2 

 Conflicts:
	src/cell.c
	src/cell.h
	src/space.c
parents fa491537 ce5d0537
Hide whitespace changes
Inline Side-by-side
examples/ExternalPointMass/run.sh
View file @ 1c1d49f4
......@@ -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
examples/IsothermalPotential/GravityOnly/README
View file @ 1c1d49f4
;
; 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.
src/cell.c
View file @ 1c1d49f4
......@@ -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);
}
src/cell.h
View file @ 1c1d49f4
......@@ -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 */
src/engine.c
View file @ 1c1d49f4
......@@ -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 */
......
src/hydro/Gizmo/hydro.h
View file @ 1c1d49f4
......@@ -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);
}
src/hydro/PressureEntropy/hydro.h
View file @ 1c1d49f4
......@@ -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);
}
......
src/hydro/PressureEntropy/hydro_iact.h
View file @ 1c1d49f4
......@@ -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;
......
src/minmax.h
View file @ 1c1d49f4
......@@ -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 */
src/riemann.h
View file @ 1c1d49f4
......@@ -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)
......
src/riemann/riemann_exact_isothermal.h 0 → 100644
View file @ 1c1d49f4
/*******************************************************************************
* 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 *