Commit 6c1f74db authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge remote-tracking branch 'origin/master' into repart-by-means

parents 08d73f87 da7892e3
......@@ -16,7 +16,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Init the project.
AC_INIT([SWIFT],[0.7.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
AC_INIT([SWIFT],[0.8.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
swift_config_flags="$*"
AC_COPYRIGHT
......
......@@ -31,7 +31,7 @@ cooling contains a temperature floor avoiding negative temperature.
Grackle
~~~~~~~
Grackle is a chemistry and cooling library presented in B. Smith et al. 2016
Grackle is a chemistry and cooling library presented in `B. Smith et al. 2016 <https://arxiv.org/abs/1610.09591>`_
(do not forget to cite if used). Four different modes are available:
equilibrium, 6 species network (H, H\\( ^+ \\), e\\( ^- \\), He, He\\( ^+ \\)
and He\\( ^{++} \\)), 9 species network (adds H\\(^-\\), H\\(_2\\) and
......
......@@ -302,6 +302,7 @@ The full section to start a DM+hydro run from Gadget DM-only ICs would
be:
.. code:: YAML
InitialConditions:
file_name: my_ics.hdf5
periodic: 1
......
......@@ -5,7 +5,7 @@ VELOCIraptor Interface
======================
This section includes information on the VELOCIraptor interface implemented in
SWIFT. There are mainly four subsection; the first section explains shortly
SWIFT. There are mainly four subsections; the first section explains shortly
how VELOCIraptor works, the second subsection explains how to configure SWIFT
with VELOCIraptor, the third subsection explains how to configure a standalone
version of VELOCIraptor and the last subsection explains how the output format
......
......@@ -46,7 +46,7 @@ The second file that is produced by VELOCIraptor is the ``.catalog_particles``
file, this file contains mainly all the IDs of the particles and has two
interesting parameters:
+ The ``Num_of_particles_in_groups`` and ``Num_of_particles_in_groups``
+ The ``Num_of_particles_in_groups`` and ``Total_num_of_particles_in_all_groups``
parameter: Gives the total number of particles in the file or the total
number of particles that are in halos.
+ The ``Particle_IDs``: The list of particles as sorted by halo, in which halo
......@@ -136,7 +136,7 @@ the unbound particles.
Properties file
---------------
The Fourth file is the ``.properties`` file, this file contains many physical
The fourth file is the ``.properties`` file, this file contains many physical
useful information of the corresponding halos. This can be divided in several
useful groups of physical parameters, on this page we have divided the several
variables which are present in the ``.properties`` file. This file has most
......@@ -179,7 +179,7 @@ Bryan and Norman 1998 properties:
+ ``Mass_BN98``, The Bryan and Norman (1998) determination of the mass of the
halo [#BN98]_.
+ ``R_BN98``, the Bryan and Norman (1998) corresponding radius[#BN98]_.
+ ``R_BN98``, the Bryan and Norman (1998) corresponding radius [#BN98]_.
Several Mass types:
"""""""""""""""""""
......
......@@ -12,7 +12,7 @@ What is VELOCIraptor?
In SWIFT it is possible to run a cosmological simulation and at the same time
do on the fly halo finding at specific predefined intervals. For finding the
Halos SWIFT uses VELOCIraptor (Elahi, Thacker and Widrow; 2011) [#velo]_, this
halos SWIFT uses VELOCIraptor (Elahi, Thacker and Widrow; 2011) [#velo]_, this
is a C++ halo finder that can use MPI. It differs from other halo finder
algorithms in the sense that it uses the velocity distributions of the
particles in the simulations and the the positions of the particles to get
......@@ -22,7 +22,7 @@ whether there are substructures in halos.
The Algorithm
-------------
The VELOCIraptor algorithm consist basically off the following steps [#ref]_:
The VELOCIraptor algorithm consists basically of the following steps [#ref]_:
1. A kd-tree is constructed based on the maximization of the Shannon-entropy,
this means that every level in the kd-tree an equal number of particles
......@@ -32,21 +32,21 @@ The VELOCIraptor algorithm consist basically off the following steps [#ref]_:
takes into account the absolute positions of the particles.
2. The next part is calculating the the centre of mass velocity and the
velocity distribution for every individual node in the kd-tree.
3. Than the algorithm estimates the background velocity density function for
3. Then the algorithm estimates the background velocity density function for
every particle based on the cell of the particle and the six nearest
neighbour cells. This prevents the background velocity density function
to be over sensitive for variations between different cells due to dominant
halo features in the velocity density function.
4. After this the algorithm searches for the nearest velocity neighbours
(:math:`N_v`) from a set of nearest position neighbours (:math:`N_x>N_v`).
The position neighbours do not need to be in the cell of the particles, in
The neighbours' positions do not need to be in the cell of the particles, in
general the set of nearest position neighbours is substantially larger than
the nearest velocity neighbours, the default is set as :math:`N_x=32 N_v`.
5. The individual local velocity density function is calculated for every
particle.
6. The fractional difference is calculated between the local velocity density
function and the background velocity density function.
7. Based on the calculated ratio outliers are picked and the outliers are
7. Based on the calculated ratio, outliers are picked and the outliers are
grouped together in halos and subhalos.
......
......@@ -7,5 +7,5 @@ then
./getIC.sh
fi
../swift -b -G -s -S -t 8 dwarf_galaxy.yml 2>&1 | tee output.log
../swift -b -G -s -S -t 8 $@ dwarf_galaxy.yml 2>&1 | tee output.log
......@@ -91,8 +91,8 @@ for i in range(402):
E_tot_snap[i] = E_kin_snap[i] + E_pot_snap[i]
Lz_snap[i] = np.sum(Lz)
print "Starting energy:", E_kin_stats[0], E_pot_stats[0], E_tot_stats[0]
print "Ending energy:", E_kin_stats[-1], E_pot_stats[-1], E_tot_stats[-1]
print("Starting energy:", E_kin_stats[0], E_pot_stats[0], E_tot_stats[0])
print("Ending energy:", E_kin_stats[-1], E_pot_stats[-1], E_tot_stats[-1])
# Plot energy evolution
figure()
......
......@@ -36,16 +36,16 @@ const_unit_length_in_cgs = (1000*PARSEC_IN_CGS)
const_unit_mass_in_cgs = (SOLAR_MASS_IN_CGS)
const_unit_velocity_in_cgs = (1e5)
print "UnitMass_in_cgs: ", const_unit_mass_in_cgs
print "UnitLength_in_cgs: ", const_unit_length_in_cgs
print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
print "UnitTime_in_cgs: ", const_unit_length_in_cgs / const_unit_velocity_in_cgs
print("UnitMass_in_cgs: ", const_unit_mass_in_cgs)
print("UnitLength_in_cgs: ", const_unit_length_in_cgs)
print("UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs)
print("UnitTime_in_cgs: ", const_unit_length_in_cgs / const_unit_velocity_in_cgs)
# derived units
const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
const_G = ((NEWTON_GRAVITY_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
print '---------------------'
print 'G in internal units: ', const_G
print('---------------------')
print('G in internal units: ', const_G)
# Parameters
......@@ -53,7 +53,7 @@ periodic = 1 # 1 For periodic box
boxSize = 100. #
max_radius = boxSize / 4. # maximum radius of particles
Mass = 1e10
print "Mass at the centre: ", Mass
print("Mass at the centre: ", Mass)
numPart = int(sys.argv[1]) # Number of particles
mass = 1.
......@@ -93,9 +93,9 @@ grp1 = file.create_group("/PartType1")
#generate particle positions
radius = max_radius * (numpy.random.rand(numPart))**(1./3.)
print '---------------------'
print 'Radius: minimum = ',min(radius)
print 'Radius: maximum = ',max(radius)
print('---------------------')
print('Radius: minimum = ',min(radius))
print('Radius: maximum = ',max(radius))
radius = numpy.sort(radius)
r = numpy.zeros((numPart, 3))
r[:,0] = radius
......@@ -104,9 +104,9 @@ r[:,0] = radius
speed = numpy.sqrt(const_G * Mass / radius)
omega = speed / radius
period = 2.*math.pi/omega
print '---------------------'
print 'Period: minimum = ',min(period)
print 'Period: maximum = ',max(period)
print('---------------------')
print('Period: minimum = ',min(period))
print('Period: maximum = ',max(period))
v = numpy.zeros((numPart, 3))
v[:,0] = -omega * r[:,1]
......
......@@ -1143,7 +1143,7 @@ int main(int argc, char *argv[]) {
/* Write final output. */
if (!force_stop) {
engine_drift_all(&e);
engine_drift_all(&e, /*drift_mpole=*/0);
engine_print_stats(&e);
#ifdef WITH_LOGGER
logger_log_all(e.logger, &e);
......
......@@ -53,8 +53,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c \
engine_marktasks.c serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
units.c common_io.c single_io.c multipole.c version.c map.c \
engine_marktasks.c engine_drift.c serial_io.c timers.c debug.c scheduler.c \
proxy.c parallel_io.c units.c common_io.c single_io.c multipole.c version.c map.c \
kernel_hydro.c tools.c part.c partition.c clocks.c parser.c \
physical_constants.c potential.c hydro_properties.c \
threadpool.c cooling.c sourceterms.c \
......
......@@ -74,6 +74,22 @@ __attribute__((always_inline)) INLINE static int cell_are_gpart_drifted(
return (c->grav.ti_old_part == e->ti_current);
}
/**
* @brief Check that the #spart in a #cell have been drifted to the current
* time.
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
* @return 1 if the #cell has been drifted to the current time, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int cell_are_spart_drifted(
const struct cell *c, const struct engine *e) {
/* Currently just use the gpart drift
* This function is just for clarity */
return cell_are_gpart_drifted(c, e);
}
/* Are cells / particles active for regular tasks ? */
/**
......@@ -185,9 +201,16 @@ __attribute__((always_inline)) INLINE static int cell_is_all_active_gravity(
__attribute__((always_inline)) INLINE static int cell_is_active_stars(
const struct cell *c, const struct engine *e) {
// LOIC: Need star-specific implementation
#ifdef SWIFT_DEBUG_CHECKS
if (c->stars.ti_end_min < e->ti_current)
error(
"cell in an impossible time-zone! c->ti_end_min=%lld (t=%e) and "
"e->ti_current=%lld (t=%e, a=%e)",
c->stars.ti_end_min, c->stars.ti_end_min * e->time_base, e->ti_current,
e->ti_current * e->time_base, e->cosmology->a);
#endif
return cell_is_active_gravity(c, e);
return (c->stars.ti_end_min == e->ti_current);
}
/**
......
......@@ -184,6 +184,7 @@ int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
pc->hydro.ti_end_max = c->hydro.ti_end_max;
pc->grav.ti_end_min = c->grav.ti_end_min;
pc->grav.ti_end_max = c->grav.ti_end_max;
pc->stars.ti_end_min = c->stars.ti_end_min;
pc->hydro.ti_old_part = c->hydro.ti_old_part;
pc->grav.ti_old_part = c->grav.ti_old_part;
pc->grav.ti_old_multipole = c->grav.ti_old_multipole;
......@@ -287,6 +288,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
c->hydro.ti_end_max = pc->hydro.ti_end_max;
c->grav.ti_end_min = pc->grav.ti_end_min;
c->grav.ti_end_max = pc->grav.ti_end_max;
c->stars.ti_end_min = pc->stars.ti_end_min;
c->hydro.ti_old_part = pc->hydro.ti_old_part;
c->grav.ti_old_part = pc->grav.ti_old_part;
c->grav.ti_old_multipole = pc->grav.ti_old_multipole;
......@@ -342,6 +344,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
temp->hydro.dx_max_part = 0.f;
temp->hydro.dx_max_sort = 0.f;
temp->stars.dx_max_part = 0.f;
temp->stars.dx_max_sort = 0.f;
temp->nodeID = c->nodeID;
temp->parent = c;
c->progeny[k] = temp;
......@@ -414,6 +417,7 @@ int cell_pack_end_step(struct cell *restrict c,
pcells[0].hydro.ti_end_max = c->hydro.ti_end_max;
pcells[0].grav.ti_end_min = c->grav.ti_end_min;
pcells[0].grav.ti_end_max = c->grav.ti_end_max;
pcells[0].stars.ti_end_min = c->stars.ti_end_min;
pcells[0].hydro.dx_max_part = c->hydro.dx_max_part;
pcells[0].stars.dx_max_part = c->stars.dx_max_part;
......@@ -451,6 +455,7 @@ int cell_unpack_end_step(struct cell *restrict c,
c->hydro.ti_end_max = pcells[0].hydro.ti_end_max;
c->grav.ti_end_min = pcells[0].grav.ti_end_min;
c->grav.ti_end_max = pcells[0].grav.ti_end_max;
c->stars.ti_end_min = pcells[0].stars.ti_end_min;
c->hydro.dx_max_part = pcells[0].hydro.dx_max_part;
c->stars.dx_max_part = pcells[0].stars.dx_max_part;
......@@ -1561,12 +1566,20 @@ void cell_check_multipole(struct cell *c) {
*/
void cell_clean(struct cell *c) {
/* Hydro */
for (int i = 0; i < 13; i++)
if (c->hydro.sort[i] != NULL) {
free(c->hydro.sort[i]);
c->hydro.sort[i] = NULL;
}
/* Stars */
for (int i = 0; i < 13; i++)
if (c->stars.sort[i] != NULL) {
free(c->stars.sort[i]);
c->stars.sort[i] = NULL;
}
/* Recurse */
for (int k = 0; k < 8; k++)
if (c->progeny[k]) cell_clean(c->progeny[k]);
......@@ -1665,13 +1678,14 @@ void cell_activate_drift_gpart(struct cell *c, struct scheduler *s) {
* @brief Activate the #spart drifts on the given cell.
*/
void cell_activate_drift_spart(struct cell *c, struct scheduler *s) {
// MATTHIEU: This will need changing
cell_activate_drift_gpart(c, s);
}
/**
* @brief Activate the sorts up a cell hierarchy.
*/
void cell_activate_sorts_up(struct cell *c, struct scheduler *s) {
void cell_activate_hydro_sorts_up(struct cell *c, struct scheduler *s) {
if (c == c->hydro.super) {
#ifdef SWIFT_DEBUG_CHECKS
......@@ -1702,7 +1716,7 @@ void cell_activate_sorts_up(struct cell *c, struct scheduler *s) {
/**
* @brief Activate the sorts on a given cell, if needed.
*/
void cell_activate_sorts(struct cell *c, int sid, struct scheduler *s) {
void cell_activate_hydro_sorts(struct cell *c, int sid, struct scheduler *s) {
/* Do we need to re-sort? */
if (c->hydro.dx_max_sort > space_maxreldx * c->dmin) {
......@@ -1711,7 +1725,7 @@ void cell_activate_sorts(struct cell *c, int sid, struct scheduler *s) {
for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
if (finger->hydro.requires_sorts) {
atomic_or(&finger->hydro.do_sort, finger->hydro.requires_sorts);
cell_activate_sorts_up(finger, s);
cell_activate_hydro_sorts_up(finger, s);
}
finger->hydro.sorted = 0;
}
......@@ -1720,7 +1734,70 @@ void cell_activate_sorts(struct cell *c, int sid, struct scheduler *s) {
/* Has this cell been sorted at all for the given sid? */
if (!(c->hydro.sorted & (1 << sid)) || c->nodeID != engine_rank) {
atomic_or(&c->hydro.do_sort, (1 << sid));
cell_activate_sorts_up(c, s);
cell_activate_hydro_sorts_up(c, s);
}
}
/**
* @brief Activate the sorts up a cell hierarchy.
*/
void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
if (c == c->super) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->stars.sorts == NULL)
error("Trying to activate un-existing c->stars.sorts");
#endif
scheduler_activate(s, c->stars.sorts);
if (c->nodeID == engine_rank) {
// MATTHIEU: to do: do we actually need both drifts here?
cell_activate_drift_part(c, s);
cell_activate_drift_spart(c, s);
}
} else {
for (struct cell *parent = c->parent;
parent != NULL && !parent->stars.do_sub_sort;
parent = parent->parent) {
parent->stars.do_sub_sort = 1;
if (parent == c->super) {
#ifdef SWIFT_DEBUG_CHECKS
if (parent->stars.sorts == NULL)
error("Trying to activate un-existing parents->stars.sorts");
#endif
scheduler_activate(s, parent->stars.sorts);
if (parent->nodeID == engine_rank) {
cell_activate_drift_part(parent, s);
cell_activate_drift_spart(parent, s);
}
break;
}
}
}
}
/**
* @brief Activate the sorts on a given cell, if needed.
*/
void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s) {
/* Do we need to re-sort? */
if (c->stars.dx_max_sort > space_maxreldx * c->dmin) {
/* Climb up the tree to active the sorts in that direction */
for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
if (finger->stars.requires_sorts) {
atomic_or(&finger->stars.do_sort, finger->stars.requires_sorts);
cell_activate_stars_sorts_up(finger, s);
}
finger->stars.sorted = 0;
}
}
/* Has this cell been sorted at all for the given sid? */
if (!(c->stars.sorted & (1 << sid)) || c->nodeID != engine_rank) {
atomic_or(&c->stars.do_sort, (1 << sid));
cell_activate_stars_sorts_up(c, s);
}
}
......@@ -2074,8 +2151,8 @@ void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
/* Do we need to sort the cells? */
cell_activate_sorts(ci, sid, s);
cell_activate_sorts(cj, sid, s);
cell_activate_hydro_sorts(ci, sid, s);
cell_activate_hydro_sorts(cj, sid, s);
}
} /* Otherwise, pair interation */
}
......@@ -2105,7 +2182,9 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
if (cj == NULL) {
/* Do anything? */
if (ci->stars.count == 0 || !cell_is_active_stars(ci, e)) return;
if (!cell_is_active_stars(ci, e) || ci->hydro.count == 0 ||
ci->stars.count == 0)
return;
/* Recurse? */
if (cell_can_recurse_in_self_stars_task(ci)) {
......@@ -2133,7 +2212,10 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
/* Should we even bother? */
if (!cell_is_active_stars(ci, e) && !cell_is_active_stars(cj, e)) return;
if (ci->stars.count == 0 || cj->stars.count == 0) return;
int should_do = ci->stars.count != 0 && cj->hydro.count != 0;
should_do |= cj->stars.count != 0 && ci->hydro.count != 0;
if (!should_do) return;
/* Get the orientation of the pair. */
double shift[3];
......@@ -2418,23 +2500,43 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
}
/* Otherwise, activate the sorts and drifts. */
else if (cell_is_active_stars(ci, e) || cell_is_active_stars(cj, e)) {
else {
/* We are going to interact this pair, so store some values. */
atomic_or(&ci->hydro.requires_sorts, 1 << sid);
atomic_or(&cj->hydro.requires_sorts, 1 << sid);
ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
if (cell_is_active_stars(ci, e) && cj->hydro.count != 0 &&
ci->stars.count != 0) {
/* We are going to interact this pair, so store some values. */
atomic_or(&cj->hydro.requires_sorts, 1 << sid);
atomic_or(&ci->stars.requires_sorts, 1 << sid);
/* Activate the drifts if the cells are local. */
if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
if (ci->nodeID == engine_rank) cell_activate_drift_spart(ci, s);
if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
if (cj->nodeID == engine_rank) cell_activate_drift_spart(cj, s);
cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
ci->stars.dx_max_sort_old = ci->stars.dx_max_sort;
/* Do we need to sort the cells? */
cell_activate_sorts(ci, sid, s);
cell_activate_sorts(cj, sid, s);
/* Activate the drifts if the cells are local. */
if (ci->nodeID == engine_rank) cell_activate_drift_spart(ci, s);
if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
/* Do we need to sort the cells? */
cell_activate_hydro_sorts(cj, sid, s);
cell_activate_stars_sorts(ci, sid, s);
}
if (cell_is_active_stars(cj, e) && ci->hydro.count != 0 &&
cj->stars.count != 0) {
/* We are going to interact this pair, so store some values. */
atomic_or(&cj->stars.requires_sorts, 1 << sid);
atomic_or(&ci->hydro.requires_sorts, 1 << sid);
ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
/* Activate the drifts if the cells are local. */
if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
if (cj->nodeID == engine_rank) cell_activate_drift_spart(cj, s);
/* Do we need to sort the cells? */
cell_activate_hydro_sorts(ci, sid, s);
cell_activate_stars_sorts(cj, sid, s);
}
}
} /* Otherwise, pair interation */
}
......@@ -2652,8 +2754,8 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
/* Check the sorts and activate them if needed. */
cell_activate_sorts(ci, t->flags, s);
cell_activate_sorts(cj, t->flags, s);
cell_activate_hydro_sorts(ci, t->flags, s);
cell_activate_hydro_sorts(cj, t->flags, s);
}
/* Store current values of dx_max and h_max. */
else if (t->type == task_type_sub_pair || t->type == task_type_sub_self) {
......@@ -2946,6 +3048,7 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
* @return 1 If the space needs rebuilding. 0 otherwise.
*/
int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
struct engine *e = s->space->e;
const int nodeID = e->nodeID;
int rebuild = 0;
......@@ -2965,25 +3068,48 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
/* Activate drifts */
if (t->type == task_type_self) {
if (ci->nodeID == nodeID) cell_activate_drift_part(ci, s);
if (ci->nodeID == nodeID) cell_activate_drift_gpart(ci, s);
if (ci->nodeID == nodeID) {
cell_activate_drift_part(ci, s);
cell_activate_drift_spart(ci, s);
}
}
/* Set the correct sorting flags and activate hydro drifts */
else if (t->type == task_type_pair) {
/* Store some values. */
atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
/* Do ci */
/* stars for ci */
atomic_or(&ci->stars.requires_sorts, 1 << t->flags);
ci->stars.dx_max_sort_old = ci->stars.dx_max_sort;
/* hydro for cj */
atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
/* Activate the drift tasks. */
if (ci->nodeID == nodeID) cell_activate_drift_part(ci, s);
if (ci->nodeID == nodeID) cell_activate_drift_spart(ci, s);
if (cj->nodeID == nodeID) cell_activate_drift_part(cj, s);
/* Check the sorts and activate them if needed. */
cell_activate_sorts(ci, t->flags, s);
cell_activate_sorts(cj, t->flags, s);
cell_activate_stars_sorts(ci, t->flags, s);
cell_activate_hydro_sorts(cj, t->flags, s);
/* Do cj */
/* hydro for ci */
atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
/* stars for cj */
atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
/* Activate the drift tasks. */
if (cj->nodeID == nodeID) cell_activate_drift_spart(cj, s);
if (ci->nodeID == nodeID) cell_activate_drift_part(ci, s);
/* Check the sorts and activate them if needed. */
cell_activate_hydro_sorts(ci, t->flags, s);
cell_activate_stars_sorts(cj, t->flags, s);
}
/* Store current values of dx_max and h_max. */
else if (t->type == task_type_sub_pair || t->type == task_type_sub_self) {
......@@ -3173,7 +3299,7 @@ void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data) {
/* All top-level cells get an MPI tag. */
#ifdef WITH_MPI
if (c->mpi.tag < 0 && c->mpi.sendto) cell_tag(c);
cell_ensure_tagged(c);
#endif
/* Super-pointer for hydro */
......@@ -3245,6 +3371,21 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
if (ti_current < ti_old_part) error("Attempt to drift to the past");
#endif
/* Early abort? */
if (c->hydro.count == 0) {
/* Clear the drift flags. */
c->hydro.do_drift = 0;
c->hydro.do_sub_drift = 0;
/* Update the time of the last drift */
c->hydro.ti_old_part = ti_current;
return;
}
/* Ok, we have some particles somewhere in the hierarchy to drift */
/* Are we not in a leaf ? */
if (c->split && (force || c->hydro.do_sub_drift)) {
......@@ -3402,6 +3543,7 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
struct spart *const sparts = c->stars.parts;