diff --git a/examples/CoolingHaloWithSpin/cooling_halo.yml b/examples/CoolingHaloWithSpin/cooling_halo.yml
index b90d2fa88146dfec89170e004c6aed5ece64bc75..684dd11fcf7adc9477d199e599dfb5b76faa91f6 100644
--- a/examples/CoolingHaloWithSpin/cooling_halo.yml
+++ b/examples/CoolingHaloWithSpin/cooling_halo.yml
@@ -10,7 +10,7 @@ InternalUnitSystem:
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 10. # The end time of the simulation (in internal units).
- dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-1 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics
@@ -32,9 +32,6 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: CoolingHalo.hdf5 # The file to read
- shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
- shift_y: 0.
- shift_z: 0.
# External potential parameters
SoftenedIsothermalPotential:
@@ -43,12 +40,12 @@ SoftenedIsothermalPotential:
position_z: 0.
vrot: 200. # rotation speed of isothermal potential in internal units
timestep_mult: 0.03 # controls time step
- epsilon: 0.1 #softening for the isothermal potential
+ epsilon: 1.0 #softening for the isothermal potential
# Cooling parameters
LambdaCooling:
- lambda_cgs: 1.0e-22 # Cooling rate (in cgs units)
+ lambda_cgs: 1.0e-22 # Cooling rate (in cgs units)
minimum_temperature: 1.0e4 # Minimal temperature (Kelvin)
mean_molecular_weight: 0.59 # Mean molecular weight
hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless)
- cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
+ cooling_tstep_mult: 0.1 # Dimensionless pre-factor for the time-step condition
diff --git a/examples/HydrostaticHalo/README b/examples/HydrostaticHalo/README
index c6630805b1e063266d4f17e9125ffaa640e5bb25..797c533926b7a131db9a08424b48b58ac3d3b8f4 100644
--- a/examples/HydrostaticHalo/README
+++ b/examples/HydrostaticHalo/README
@@ -1,10 +1,12 @@
+Hydrostatic halo in equilibrium in an isothermal potential. Running
+for 10 dynamical times.
To make the initial conditions we distribute gas particles randomly in
a cube with a side length twice that of the virial radius. The density
profile of the gas is proportional to r^(-2) where r is the distance
from the centre of the cube.
-The parameter v_rot (in makeIC.py and cooling.yml) sets the circular
+The parameter v_rot (in makeIC.py and hydrostatic.yml) sets the circular
velocity of the halo, and by extension, the viral radius, viral mass,
and the internal energy of the gas such that hydrostatic equilibrium
is achieved.
@@ -12,10 +14,12 @@ is achieved.
To run this example, make such that the code is compiled with either
the isothermal potential or softened isothermal potential set in
src/const.h. In the latter case, a (small) value of epsilon needs to
-be set in cooling.yml. 0.1 kpc should work well.
+be set in hydrostatic.yml. ~1 kpc should work well.
The plotting scripts produce a plot of the density, internal energy
-and radial velocity profile for each
-snapshot. test_energy_conservation.py shows the evolution of energy
+and radial velocity profile for each snapshot and divides the profile
+by the expected profile.
+
+The script test_energy_conservation.py shows the evolution of energy
with time. These can be used to check if the example has run properly.
diff --git a/examples/HydrostaticHalo/hydrostatic.yml b/examples/HydrostaticHalo/hydrostatic.yml
index f0f07f943c768181cda459f79dc31bf18ca8d20a..39a91a4ec475a70ef4e61b9cdc59b8221a74093e 100644
--- a/examples/HydrostaticHalo/hydrostatic.yml
+++ b/examples/HydrostaticHalo/hydrostatic.yml
@@ -9,7 +9,7 @@ InternalUnitSystem:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1.0 # The end time of the simulation (in internal units).
+ time_end: 30. # The end time of the simulation (in internal units).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
@@ -20,21 +20,18 @@ Statistics:
# Parameters governing the snapshots
Snapshots:
basename: Hydrostatic # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
- delta_neighbours: 1. # The tolerance for the targetted number of neighbours.
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions
InitialConditions:
file_name: Hydrostatic.hdf5 # The file to read
- shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
- shift_y: 0.
- shift_z: 0.
# External potential parameters
SoftenedIsothermalPotential:
@@ -42,6 +39,6 @@ SoftenedIsothermalPotential:
position_y: 0.
position_z: 0.
vrot: 200. # rotation speed of isothermal potential in internal units
- epsilon: 0.1
+ epsilon: 1.0
timestep_mult: 0.03 # controls time step
diff --git a/examples/HydrostaticHalo/makeIC.py b/examples/HydrostaticHalo/makeIC.py
index 376423fa65131761594cfa25c3beb3252c3750fb..f33387e18dd0ab523684227f5c745b5c8b807b7f 100644
--- a/examples/HydrostaticHalo/makeIC.py
+++ b/examples/HydrostaticHalo/makeIC.py
@@ -116,9 +116,9 @@ print "x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0]))
print "y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1]))
print "z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2]))
-print np.mean(coords[:,0])
-print np.mean(coords[:,1])
-print np.mean(coords[:,2])
+#print np.mean(coords[:,0])
+#print np.mean(coords[:,1])
+#print np.mean(coords[:,2])
#now find the particles which are within the box
diff --git a/examples/HydrostaticHalo/run.sh b/examples/HydrostaticHalo/run.sh
index f3d3694760ced9e41719de64229f0b5565adf5d4..d23ead6a67f43c9d19d76a797e72d050a3978d61 100755
--- a/examples/HydrostaticHalo/run.sh
+++ b/examples/HydrostaticHalo/run.sh
@@ -2,14 +2,23 @@
# Generate the initial conditions if they are not present.
echo "Generating initial conditions for the isothermal potential box example..."
-python makeIC.py 10000
+python makeIC.py 100000
-../swift -g -s -t 16 hydrostatic.yml 2>&1 | tee output.log
+# Run for 10 dynamical times
+../swift -g -s -t 2 hydrostatic.yml 2>&1 | tee output.log
-python density_profile.py 2. 200 100
+echo "Plotting density profiles"
+mkdir plots
+mkdir plots/density_profile
+python density_profile.py 2. 200 300
-python internal_energy_profile.py 2. 200 100
+echo "Plotting internal energy profiles"
+mkdir plots/internal_energy
+python internal_energy_profile.py 2. 200 300
-python velocity_profile.py 2. 200 100
+echo "Plotting radial velocity profiles"
+mkdir plots/radial_velocity_profile
+python velocity_profile.py 2. 200 300
-python test_energy_conservation.py
+echo "Plotting energy as a function of time"
+python test_energy_conservation.py 300
diff --git a/examples/HydrostaticHalo/test_energy_conservation.py b/examples/HydrostaticHalo/test_energy_conservation.py
index 26041890c0e5c8b80abda4448f0ed55d8e3b98d4..ca091050c4127d11a37a2cc7504e42d244031e25 100644
--- a/examples/HydrostaticHalo/test_energy_conservation.py
+++ b/examples/HydrostaticHalo/test_energy_conservation.py
@@ -3,7 +3,7 @@ import h5py as h5
import matplotlib.pyplot as plt
import sys
-n_snaps = 5
+n_snaps = int(sys.argv[1])
#some constants
OMEGA = 0.3 # Cosmological matter fraction at z = 0
@@ -24,7 +24,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["IsothermalPotential:vrot"])
+v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
v_c_cgs = v_c * unit_velocity_cgs
header = f["Header"]
N = header.attrs["NumPart_Total"][0]
diff --git a/examples/IsothermalPotential/GravityOnly/README b/examples/IsothermalPotential/README
similarity index 100%
rename from examples/IsothermalPotential/GravityOnly/README
rename to examples/IsothermalPotential/README
diff --git a/examples/IsothermalPotential/GravityOnly/isothermal.yml b/examples/IsothermalPotential/isothermal.yml
similarity index 100%
rename from examples/IsothermalPotential/GravityOnly/isothermal.yml
rename to examples/IsothermalPotential/isothermal.yml
diff --git a/examples/IsothermalPotential/GravityOnly/makeIC.py b/examples/IsothermalPotential/makeIC.py
similarity index 97%
rename from examples/IsothermalPotential/GravityOnly/makeIC.py
rename to examples/IsothermalPotential/makeIC.py
index 07993f19d40a9a3b9a4b86c9dd8c44f7e6fa3d7e..976119f0a312c5acc81fab943ba3cf5769102269 100644
--- a/examples/IsothermalPotential/GravityOnly/makeIC.py
+++ b/examples/IsothermalPotential/makeIC.py
@@ -141,17 +141,17 @@ ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
v = numpy.zeros(1)
-m = numpy.full((numPart, ), mass)
+m = numpy.full((numPart, ), mass, dtype='f')
ds = grp1.create_dataset('Masses', (numPart,), 'f')
ds[()] = m
m = numpy.zeros(1)
-h = numpy.full((numPart, ), 1.1255 * boxSize / L)
+h = numpy.full((numPart, ), 1.1255 * boxSize / L, dtype='f')
ds = grp1.create_dataset('SmoothingLength', (numPart,), 'f')
ds[()] = h
h = numpy.zeros(1)
-u = numpy.full((numPart, ), internalEnergy)
+u = numpy.full((numPart, ), internalEnergy, dtype='f')
ds = grp1.create_dataset('InternalEnergy', (numPart,), 'f')
ds[()] = u
u = numpy.zeros(1)
diff --git a/examples/IsothermalPotential/GravityOnly/run.sh b/examples/IsothermalPotential/run.sh
similarity index 79%
rename from examples/IsothermalPotential/GravityOnly/run.sh
rename to examples/IsothermalPotential/run.sh
index f6adfcececf4923485c0deabd97e9af9a6f64b05..28a3cc0910f986f84bcd603091543643356f1c4a 100755
--- a/examples/IsothermalPotential/GravityOnly/run.sh
+++ b/examples/IsothermalPotential/run.sh
@@ -7,4 +7,4 @@ then
python makeIC.py 1000 1
fi
-../../swift -g -t 2 isothermal.yml 2>&1 | tee output.log
+../swift -g -t 1 isothermal.yml 2>&1 | tee output.log
diff --git a/examples/IsothermalPotential/GravityOnly/test.pro b/examples/IsothermalPotential/test.pro
similarity index 100%
rename from examples/IsothermalPotential/GravityOnly/test.pro
rename to examples/IsothermalPotential/test.pro
diff --git a/examples/Makefile.am b/examples/Makefile.am
index f0c37bed62e6c95c6210cf48c55dae2295e338de..4da84788a485dacd2103fe85ad3e729ade6b582a 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -65,6 +65,9 @@ EXTRA_DIST = BigCosmoVolume/makeIC.py \
EAGLE_50/eagle_50.yml EAGLE_50/getIC.sh EAGLE_50/README EAGLE_50/run.sh \
ExternalPointMass/externalPointMass.yml ExternalPointMass/makeIC.py ExternalPointMass/run.sh ExternalPointMass/test.pro \
GreshoVortex_2D/getGlass.sh GreshoVortex_2D/gresho.yml GreshoVortex_2D/makeIC.py GreshoVortex_2D/plotSolution.py GreshoVortex_2D/run.sh \
+ HydrostaticHalo/README HydrostaticHalo/hydrostatic.yml HydrostaticHalo/makeIC.py HydrostaticHalo/run.sh \
+ HydrostaticHalo/density_profile.py HydrostaticHalo/velocity_profile.py HydrostaticHalo/internal_energy_profile.py HydrostaticHalo/test_energy_conservation.py \
+ IsothermalPotential/README IsothermalPotential/run.sh IsothermalPotential/test.pro IsothermalPotential/isothermal.yml IsothermalPotential/makeIC.py \
KelvinHelmholtz_2D/kelvinHelmholtz.yml KelvinHelmholtz_2D/makeIC.py KelvinHelmholtz_2D/plotSolution.py KelvinHelmholtz_2D/run.sh \
MultiTypes/makeIC.py MultiTypes/multiTypes.yml MultiTypes/run.sh \
PerturbedBox_2D/makeIC.py PerturbedBox_2D/perturbedPlane.yml \
diff --git a/m4/ax_cc_maxopt.m4 b/m4/ax_cc_maxopt.m4
index d5b87c903c5af46b767d7b49e40963345c8128af..93d5d6dcd78ff77c934f77ad0e1e02ef37873a37 100644
--- a/m4/ax_cc_maxopt.m4
+++ b/m4/ax_cc_maxopt.m4
@@ -160,6 +160,9 @@ if test "$ac_test_CFLAGS" != "set"; then
# note that we enable "unsafe" fp optimization with other compilers, too
AX_CHECK_COMPILE_FLAG(-ffast-math, CFLAGS="$CFLAGS -ffast-math")
+ # not all codes will benefit from this.
+ AX_CHECK_COMPILE_FLAG(-funroll-loops, CFLAGS="$CFLAGS -funroll-loops")
+
AX_GCC_ARCHFLAG($acx_maxopt_portable)
;;
diff --git a/src/Makefile.am b/src/Makefile.am
index 30afe6cf16960fe5b13f7a0bf6cb0cae86eaa43b..8b8dce69877b1c561d96b6d5720fdab1d068436b 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -60,7 +60,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \
units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
- dimension.h equation_of_state.h \
+ dimension.h equation_of_state.h active.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
diff --git a/src/active.h b/src/active.h
new file mode 100644
index 0000000000000000000000000000000000000000..17adfd07d2bfe0519f0b432217d5253de13c4b78
--- /dev/null
+++ b/src/active.h
@@ -0,0 +1,104 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_ACTIVE_H
+#define SWIFT_ACTIVE_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes. */
+#include "cell.h"
+#include "const.h"
+#include "engine.h"
+#include "part.h"
+
+/**
+ * @brief Does a cell contain any active particle ?
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ */
+__attribute__((always_inline)) INLINE static int cell_is_active(
+ const struct cell *c, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->ti_end_min < e->ti_current)
+ error("cell in an impossible time-zone! c->ti_end_min=%d e->ti_current=%d",
+ c->ti_end_min, e->ti_current);
+#endif
+
+ return (c->ti_end_min == e->ti_current);
+}
+
+/**
+ * @brief Are *all* particles in a cell active ?
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ */
+__attribute__((always_inline)) INLINE static int cell_is_all_active(
+ const struct cell *c, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->ti_end_max < e->ti_current)
+ error("cell in an impossible time-zone! c->ti_end_max=%d e->ti_current=%d",
+ c->ti_end_max, e->ti_current);
+#endif
+
+ return (c->ti_end_max == e->ti_current);
+}
+
+/**
+ * @brief Is this particle active ?
+ *
+ * @param p The #part.
+ * @param e The #engine containing information about the current time.
+ */
+__attribute__((always_inline)) INLINE static int part_is_active(
+ const struct part *p, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (p->ti_end < e->ti_current)
+ error("particle in an impossible time-zone! p->ti_end=%d e->ti_current=%d",
+ p->ti_end, e->ti_current);
+#endif
+
+ return (p->ti_end == e->ti_current);
+}
+
+/**
+ * @brief Is this g-particle active ?
+ *
+ * @param gp The #gpart.
+ * @param e The #engine containing information about the current time.
+ */
+__attribute__((always_inline)) INLINE static int gpart_is_active(
+ const struct gpart *gp, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (gp->ti_end < e->ti_current)
+ error(
+ "g-particle in an impossible time-zone! gp->ti_end=%d e->ti_current=%d",
+ gp->ti_end, e->ti_current);
+#endif
+
+ return (gp->ti_end == e->ti_current);
+}
+
+#endif /* SWIFT_ACTIVE_H */
diff --git a/src/cell.c b/src/cell.c
index 573272d05839d6d082dac61c97f6abd18d8eb41a..18c2279b76cdd54a626dfe91f655c9e4e01d6b1b 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -47,6 +47,7 @@
#include "cell.h"
/* Local headers. */
+#include "active.h"
#include "atomic.h"
#include "error.h"
#include "gravity.h"
@@ -867,14 +868,14 @@ void cell_clean(struct cell *c) {
* @brief Checks whether a given cell needs drifting or not.
*
* @param c the #cell.
- * @param ti_current The current time on the integer time-line.
+ * @param e The #engine (holding current time information).
*
* @return 1 If the cell needs drifting, 0 otherwise.
*/
-int cell_is_drift_needed(struct cell *c, int ti_current) {
+int cell_is_drift_needed(struct cell *c, const struct engine *e) {
/* Do we have at least one active particle in the cell ?*/
- if (c->ti_end_min == ti_current) return 1;
+ if (cell_is_active(c, e)) return 1;
/* Loop over the pair tasks that involve this cell */
for (struct link *l = c->density; l != NULL; l = l->next) {
@@ -882,9 +883,9 @@ int cell_is_drift_needed(struct cell *c, int ti_current) {
if (l->t->type != task_type_pair && l->t->type != task_type_sub_pair)
continue;
- /* Does the other cell in the pair have an active particle ? */
- if ((l->t->ci == c && l->t->cj->ti_end_min == ti_current) ||
- (l->t->cj == c && l->t->ci->ti_end_min == ti_current))
+ /* Is the other cell in the pair active ? */
+ if ((l->t->ci == c && cell_is_active(l->t->cj, e)) ||
+ (l->t->cj == c && cell_is_active(l->t->ci, e)))
return 1;
}
diff --git a/src/cell.h b/src/cell.h
index 9e5bed091178b59e1b757c420bc1d5fde0b9ce42..10c67be768352374b47736bed225a136f05dba28 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -37,6 +37,7 @@
#include "task.h"
/* Avoid cyclic inclusions */
+struct engine;
struct space;
struct scheduler;
@@ -290,7 +291,7 @@ int cell_are_neighbours(const struct cell *restrict ci,
const struct cell *restrict cj);
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_is_drift_needed(struct cell *c, const struct engine *e);
int cell_unskip_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index 0636adec588dd2d215f9aaf2874af61a707ab532..cb9db2dc34a6014ea15a24d368a006fee3838d67 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -24,6 +24,7 @@
#define SWIFT_COOLING_CONST_LAMBDA_H
/* Some standard headers. */
+#include <float.h>
#include <math.h>
/* Local includes. */
@@ -84,7 +85,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
/* Calculate du_dt */
const float du_dt = cooling_rate(phys_const, us, cooling, p);
- /* Intergrate cooling equation, but enforce energy floor */
+ /* Integrate cooling equation, but enforce energy floor */
float u_new;
if (u_old + du_dt * dt > u_floor) {
u_new = u_old + du_dt * dt;
@@ -92,6 +93,9 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
u_new = u_floor;
}
+ /* Don't allow particle to cool too much in one timestep */
+ if (u_new < 0.5f * u_old) u_new = 0.5f * u_old;
+
/* Update the internal energy */
hydro_set_internal_energy(p, u_new);
@@ -112,13 +116,16 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
const struct phys_const* restrict phys_const,
const struct UnitSystem* restrict us, const struct part* restrict p) {
- /* Get du_dt */
- const float du_dt = cooling_rate(phys_const, us, cooling, p);
-
/* Get current internal energy (dt=0) */
const float u = hydro_get_internal_energy(p, 0.f);
+ const float du_dt = cooling_rate(phys_const, us, cooling, p);
- return u / fabsf(du_dt);
+ /* If we are close to (or below) the energy floor, we ignore cooling timestep
+ */
+ if (u < 1.01f * cooling->min_energy)
+ return FLT_MAX;
+ else
+ return cooling->cooling_tstep_mult * u / fabsf(du_dt);
}
/**
diff --git a/src/debug.c b/src/debug.c
index 196c2f7e49afa827f8d80539cdeb3a975cbf31fc..5be0370f64f2c21e6b28c40cc9802520087ae07f 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -193,6 +193,59 @@ int checkSpacehmax(struct space *s) {
return 0;
}
+/**
+ * @brief Check if the h_max and dx_max values of a cell's hierarchy are
+ * consistent with the particles. Report verbosely if not.
+ *
+ * @param c the top cell of the hierarchy.
+ * @param depth the recursion depth for use in messages. Set to 0 initially.
+ * @result 1 or 0
+ */
+int checkCellhdxmax(const struct cell *c, int *depth) {
+
+ *depth = *depth + 1;
+
+ float h_max = 0.0f;
+ float dx_max = 0.0f;
+ if (!c->split) {
+ const size_t nr_parts = c->count;
+ struct part *parts = c->parts;
+ for (size_t k = 0; k < nr_parts; k++) {
+ h_max = (h_max > parts[k].h) ? h_max : parts[k].h;
+ }
+ } else {
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) {
+ struct cell *cp = c->progeny[k];
+ checkCellhdxmax(cp, depth);
+ dx_max = max(dx_max, cp->dx_max);
+ h_max = max(h_max, cp->h_max);
+ }
+ }
+
+ /* Check. */
+ int result = 1;
+ if (c->h_max != h_max) {
+ message("%d Inconsistent h_max: cell %f != parts %f", *depth, c->h_max, h_max);
+ message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+ result = 0;
+ }
+ if (c->dx_max != dx_max) {
+ message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max, dx_max);
+ message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+ result = 0;
+ }
+
+ /* Check rebuild criterion. */
+ if (h_max > c->dmin) {
+ message("%d Inconsistent c->dmin: %f > %f", *depth, h_max, c->dmin);
+ message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+ result = 0;
+ }
+
+ return result;
+}
+
#ifdef HAVE_METIS
/**
diff --git a/src/debug.h b/src/debug.h
index d75ce91f1ea23221626f32e11472f713e9731789..7422a6f7f9815490966f08415e0312876ce0123f 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -32,6 +32,7 @@ void printParticle_single(const struct part *p, const struct xpart *xp);
void printgParticle_single(struct gpart *gp);
int checkSpacehmax(struct space *s);
+int checkCellhdxmax(const struct cell *c, int *depth);
#ifdef HAVE_METIS
#include "metis.h"
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index ccdd0cee32b9386eff54da655b75285b8e08a598..3fd357a2d8778f5ca8b014935d538350eccb99c6 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -199,10 +199,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* and add the self-contribution term.
*
* @param p The particle to act upon
- * @param time The current time
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
- struct part *restrict p, float time) {
+ struct part *restrict p) {
/* Some smoothing length multiples. */
const float h = p->h;
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index c999e20d401570ac6518291df8cf315569fe78bd..157893bc9e27806d2b97ac5f5a81d0f6fbb1c589 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -246,10 +246,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* and add the self-contribution term.
*
* @param p The particle to act upon
- * @param ti_current The current time (on the integer timeline)
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
- struct part *restrict p, int ti_current) {
+ struct part *restrict p) {
/* Some smoothing length multiples. */
const float h = p->h;
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index bd970795bdf070a7bd7915cc4f493218dbf319d1..1c64291ee64dd770b1f1a76371f67a34230365c7 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -109,10 +109,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* passive particles.
*
* @param p The particle to act upon.
- * @param The current physical time.
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
- struct part* restrict p, float time) {
+ struct part* restrict p) {
/* Some smoothing length multiples. */
const float h = p->h;
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 3015f26c6bad7f006e8cda16695c750ff40d74df..3b3454f1bb348b178ac57899da4f7611802a69cd 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -256,10 +256,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* added to them here.
*
* @param p The particle to act upon
- * @param time The current time
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
- struct part *restrict p, float time) {
+ struct part *restrict p) {
/* Some smoothing length multiples. */
const float h = p->h;
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 0c69728a9ce6317a8d7ddab945e7aab6e94ee247..8c063596efd3be97ebb4da6b6879ac06122bd357 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -254,10 +254,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* and add the self-contribution term.
*
* @param p The particle to act upon
- * @param ti_current The current time (on the integer timeline)
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
- struct part *restrict p, int ti_current) {
+ struct part *restrict p) {
/* Some smoothing length multiples. */
const float h = p->h;
diff --git a/src/runner.c b/src/runner.c
index 472f11c067d3e83018d58d200cff2ca0271eb3d2..2d6da4e4aedc9c40d1dade243e605e9aeda86dbe 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -38,6 +38,7 @@
#include "runner.h"
/* Local headers. */
+#include "active.h"
#include "approx_math.h"
#include "atomic.h"
#include "cell.h"
@@ -132,16 +133,16 @@ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_sourceterms(r, c->progeny[k], 0);
- return;
- }
+ } else {
- if (count > 0) {
+ if (count > 0) {
- /* do sourceterms in this cell? */
- const int incell =
- sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
- if (incell == 1) {
- sourceterms_apply(r, sourceterms, c);
+ /* do sourceterms in this cell? */
+ const int incell =
+ sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
+ if (incell == 1) {
+ sourceterms_apply(r, sourceterms, c);
+ }
}
}
@@ -159,37 +160,32 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gparts = c->gparts;
const int gcount = c->gcount;
- const int ti_current = r->e->ti_current;
- const struct external_potential *potential = r->e->external_potential;
- const struct phys_const *constants = r->e->physical_constants;
+ const struct engine *e = r->e;
+ const struct external_potential *potential = e->external_potential;
+ const struct phys_const *constants = e->physical_constants;
const double time = r->e->time;
TIMER_TIC;
/* Anything to do here? */
- if (c->ti_end_min > ti_current) return;
+ if (!cell_is_active(c, e)) return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_grav_external(r, c->progeny[k], 0);
- return;
- }
-
-#ifdef TASK_VERBOSE
- OUT;
-#endif
-
- /* Loop over the gparts in this cell. */
- for (int i = 0; i < gcount; i++) {
+ } else {
- /* Get a direct pointer on the part. */
- struct gpart *restrict gp = &gparts[i];
+ /* Loop over the gparts in this cell. */
+ for (int i = 0; i < gcount; i++) {
- /* Is this part within the time step? */
- if (gp->ti_end <= ti_current) {
+ /* Get a direct pointer on the part. */
+ struct gpart *restrict gp = &gparts[i];
- external_gravity_acceleration(time, potential, constants, gp);
+ /* Is this part within the time step? */
+ if (gpart_is_active(gp, e)) {
+ external_gravity_acceleration(time, potential, constants, gp);
+ }
}
}
@@ -221,26 +217,22 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
- return;
- }
-
-#ifdef TASK_VERBOSE
- OUT;
-#endif
+ } else {
- /* Loop over the parts in this cell. */
- for (int i = 0; i < count; i++) {
+ /* Loop over the parts in this cell. */
+ for (int i = 0; i < count; i++) {
- /* Get a direct pointer on the part. */
- struct part *restrict p = &parts[i];
- struct xpart *restrict xp = &xparts[i];
+ /* Get a direct pointer on the part. */
+ struct part *restrict p = &parts[i];
+ struct xpart *restrict xp = &xparts[i];
- /* Kick has already updated ti_end, so need to check ti_begin */
- if (p->ti_begin == ti_current) {
+ /* Kick has already updated ti_end, so need to check ti_begin */
+ if (p->ti_begin == ti_current) {
- const double dt = (p->ti_end - p->ti_begin) * timeBase;
+ const double dt = (p->ti_end - p->ti_begin) * timeBase;
- cooling_cool_part(constants, us, cooling_func, p, xp, dt);
+ cooling_cool_part(constants, us, cooling_func, p, xp, dt);
+ }
}
}
@@ -492,18 +484,17 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gparts = c->gparts;
const int count = c->count;
const int gcount = c->gcount;
- const int ti_current = r->e->ti_current;
+ const struct engine *e = r->e;
TIMER_TIC;
/* Anything to do here? */
- if (c->ti_end_min > ti_current) return;
+ if (!cell_is_active(c, e)) return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_init(r, c->progeny[k], 0);
- return;
} else {
/* Loop over the parts in this cell. */
@@ -512,7 +503,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
/* Get a direct pointer on the part. */
struct part *restrict p = &parts[i];
- if (p->ti_end <= ti_current) {
+ if (part_is_active(p, e)) {
/* Get ready for a density calculation */
hydro_init_part(p);
@@ -525,7 +516,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
/* Get a direct pointer on the part. */
struct gpart *restrict gp = &gparts[i];
- if (gp->ti_end <= ti_current) {
+ if (gpart_is_active(gp, e)) {
/* Get ready for a density calculation */
gravity_init_gpart(gp);
@@ -542,23 +533,25 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
*
* @param r The runner thread.
* @param c The cell.
+ * @param timer Are we timing this ?
*/
-void runner_do_extra_ghost(struct runner *r, struct cell *c) {
+void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
#ifdef EXTRA_HYDRO_LOOP
struct part *restrict parts = c->parts;
const int count = c->count;
- const int ti_current = r->e->ti_current;
+ const struct engine *e = r->e;
+
+ TIMER_TIC;
/* Anything to do here? */
- if (c->ti_end_min > ti_current) return;
+ if (!cell_is_active(c, e)) return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k]);
- return;
+ if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k], 0);
} else {
/* Loop over the parts in this cell. */
@@ -567,7 +560,7 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c) {
/* Get a direct pointer on the part. */
struct part *restrict p = &parts[i];
- if (p->ti_end <= ti_current) {
+ if (part_is_active(p, e)) {
/* Get ready for a force calculation */
hydro_end_gradient(p);
@@ -575,6 +568,8 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c) {
}
}
+ if (timer) TIMER_TOC(timer_do_extra_ghost);
+
#else
error("SWIFT was not compiled with the extra hydro loop activated.");
#endif
@@ -586,164 +581,166 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c) {
*
* @param r The runner thread.
* @param c The cell.
+ * @param timer Are we timing this ?
*/
-void runner_do_ghost(struct runner *r, struct cell *c) {
+void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
int redo, count = c->count;
- const int ti_current = r->e->ti_current;
- const double timeBase = r->e->timeBase;
- const float target_wcount = r->e->hydro_properties->target_neighbours;
+ const struct engine *e = r->e;
+ const int ti_current = e->ti_current;
+ const double timeBase = e->timeBase;
+ const float target_wcount = e->hydro_properties->target_neighbours;
const float max_wcount =
- target_wcount + r->e->hydro_properties->delta_neighbours;
+ target_wcount + e->hydro_properties->delta_neighbours;
const float min_wcount =
- target_wcount - r->e->hydro_properties->delta_neighbours;
- const int max_smoothing_iter =
- r->e->hydro_properties->max_smoothing_iterations;
+ target_wcount - e->hydro_properties->delta_neighbours;
+ const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
TIMER_TIC;
/* Anything to do here? */
- if (c->ti_end_min > ti_current) return;
+ if (!cell_is_active(c, e)) return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k]);
- return;
- }
+ if (c->progeny[k] != NULL) runner_do_ghost(r, c->progeny[k], 0);
+ } else {
- /* Init the IDs that have to be updated. */
- int *pid = NULL;
- if ((pid = malloc(sizeof(int) * count)) == NULL)
- error("Can't allocate memory for pid.");
- for (int k = 0; k < count; k++) pid[k] = k;
+ /* Init the IDs that have to be updated. */
+ int *pid = NULL;
+ if ((pid = malloc(sizeof(int) * count)) == NULL)
+ error("Can't allocate memory for pid.");
+ for (int k = 0; k < count; k++) pid[k] = k;
- /* While there are particles that need to be updated... */
- for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
- num_reruns++) {
+ /* While there are particles that need to be updated... */
+ for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
+ num_reruns++) {
- /* Reset the redo-count. */
- redo = 0;
+ /* Reset the redo-count. */
+ redo = 0;
- /* Loop over the parts in this cell. */
- for (int i = 0; i < count; i++) {
+ /* Loop over the parts in this cell. */
+ for (int i = 0; i < count; i++) {
- /* Get a direct pointer on the part. */
- struct part *restrict p = &parts[pid[i]];
- struct xpart *restrict xp = &xparts[pid[i]];
+ /* Get a direct pointer on the part. */
+ struct part *restrict p = &parts[pid[i]];
+ struct xpart *restrict xp = &xparts[pid[i]];
- /* Is this part within the timestep? */
- if (p->ti_end <= ti_current) {
+ /* Is this part within the timestep? */
+ if (part_is_active(p, e)) {
- /* Finish the density calculation */
- hydro_end_density(p, ti_current);
+ /* Finish the density calculation */
+ hydro_end_density(p);
- float h_corr = 0.f;
+ float h_corr = 0.f;
- /* If no derivative, double the smoothing length. */
- if (p->density.wcount_dh == 0.0f) h_corr = p->h;
+ /* If no derivative, double the smoothing length. */
+ if (p->density.wcount_dh == 0.0f) h_corr = p->h;
- /* Otherwise, compute the smoothing length update (Newton step). */
- else {
- h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
+ /* Otherwise, compute the smoothing length update (Newton step). */
+ else {
+ h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
- /* Truncate to the range [ -p->h/2 , p->h ]. */
- h_corr = (h_corr < p->h) ? h_corr : p->h;
- h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
- }
+ /* Truncate to the range [ -p->h/2 , p->h ]. */
+ h_corr = (h_corr < p->h) ? h_corr : p->h;
+ h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
+ }
- /* Did we get the right number density? */
- if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
+ /* Did we get the right number density? */
+ if (p->density.wcount > max_wcount ||
+ p->density.wcount < min_wcount) {
- /* Ok, correct then */
- p->h += h_corr;
+ /* Ok, correct then */
+ p->h += h_corr;
- /* Flag for another round of fun */
- pid[redo] = pid[i];
- redo += 1;
+ /* Flag for another round of fun */
+ pid[redo] = pid[i];
+ redo += 1;
- /* Re-initialise everything */
- hydro_init_part(p);
+ /* Re-initialise everything */
+ hydro_init_part(p);
- /* Off we go ! */
- continue;
- }
+ /* Off we go ! */
+ continue;
+ }
- /* We now have a particle whose smoothing length has converged */
+ /* We now have a particle whose smoothing length has converged */
- /* As of here, particle force variables will be set. */
+ /* As of here, particle force variables will be set. */
- /* Compute variables required for the force loop */
- hydro_prepare_force(p, xp, ti_current, timeBase);
+ /* Compute variables required for the force loop */
+ hydro_prepare_force(p, xp, ti_current, timeBase);
- /* The particle force values are now set. Do _NOT_
- try to read any particle density variables! */
+ /* The particle force values are now set. Do _NOT_
+ try to read any particle density variables! */
- /* Prepare the particle for the force loop over neighbours */
- hydro_reset_acceleration(p);
+ /* Prepare the particle for the force loop over neighbours */
+ hydro_reset_acceleration(p);
+ }
}
- }
- /* We now need to treat the particles whose smoothing length had not
- * converged again */
+ /* We now need to treat the particles whose smoothing length had not
+ * converged again */
- /* Re-set the counter for the next loop (potentially). */
- count = redo;
- if (count > 0) {
+ /* Re-set the counter for the next loop (potentially). */
+ count = redo;
+ if (count > 0) {
- /* Climb up the cell hierarchy. */
- for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
+ /* Climb up the cell hierarchy. */
+ for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
- /* Run through this cell's density interactions. */
- for (struct link *l = finger->density; l != NULL; l = l->next) {
+ /* Run through this cell's density interactions. */
+ for (struct link *l = finger->density; l != NULL; l = l->next) {
- /* Self-interaction? */
- if (l->t->type == task_type_self)
- runner_doself_subset_density(r, finger, parts, pid, count);
+ /* Self-interaction? */
+ if (l->t->type == task_type_self)
+ runner_doself_subset_density(r, finger, parts, pid, count);
- /* Otherwise, pair interaction? */
- else if (l->t->type == task_type_pair) {
+ /* Otherwise, pair interaction? */
+ else if (l->t->type == task_type_pair) {
- /* Left or right? */
- if (l->t->ci == finger)
- runner_dopair_subset_density(r, finger, parts, pid, count,
- l->t->cj);
- else
- runner_dopair_subset_density(r, finger, parts, pid, count,
- l->t->ci);
+ /* Left or right? */
+ if (l->t->ci == finger)
+ runner_dopair_subset_density(r, finger, parts, pid, count,
+ l->t->cj);
+ else
+ runner_dopair_subset_density(r, finger, parts, pid, count,
+ l->t->ci);
- }
+ }
- /* Otherwise, sub-self interaction? */
- else if (l->t->type == task_type_sub_self)
- runner_dosub_subset_density(r, finger, parts, pid, count, NULL, -1,
- 1);
-
- /* Otherwise, sub-pair interaction? */
- else if (l->t->type == task_type_sub_pair) {
-
- /* Left or right? */
- if (l->t->ci == finger)
- runner_dosub_subset_density(r, finger, parts, pid, count,
- l->t->cj, -1, 1);
- else
- runner_dosub_subset_density(r, finger, parts, pid, count,
- l->t->ci, -1, 1);
+ /* Otherwise, sub-self interaction? */
+ else if (l->t->type == task_type_sub_self)
+ runner_dosub_subset_density(r, finger, parts, pid, count, NULL,
+ -1, 1);
+
+ /* Otherwise, sub-pair interaction? */
+ else if (l->t->type == task_type_sub_pair) {
+
+ /* Left or right? */
+ if (l->t->ci == finger)
+ runner_dosub_subset_density(r, finger, parts, pid, count,
+ l->t->cj, -1, 1);
+ else
+ runner_dosub_subset_density(r, finger, parts, pid, count,
+ l->t->ci, -1, 1);
+ }
}
}
}
}
- }
- if (count)
- message("Smoothing length failed to converge on %i particles.", count);
+ if (count)
+ message("Smoothing length failed to converge on %i particles.", count);
- /* Be clean */
- free(pid);
+ /* Be clean */
+ free(pid);
+ }
- TIMER_TOC(timer_do_ghost);
+ if (timer) TIMER_TOC(timer_do_ghost);
}
/**
@@ -757,20 +754,18 @@ void runner_do_ghost(struct runner *r, struct cell *c) {
*/
static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
- const int ti_current = e->ti_current;
-
/* Unskip any active tasks. */
- if (c->ti_end_min == e->ti_current) {
+ if (cell_is_active(c, e)) {
const int forcerebuild = cell_unskip_tasks(c, &e->sched);
if (forcerebuild) atomic_inc(&e->forcerebuild);
}
/* Do we really need to drift? */
if (drift) {
- if (!e->drift_all && !cell_is_drift_needed(c, ti_current)) return;
+ if (!e->drift_all && !cell_is_drift_needed(c, e)) return;
} else {
- /* Not drifting, but may still need to recurse for task skipping. */
+ /* Not drifting, but may still need to recurse for task un-skipping. */
if (c->split) {
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
@@ -782,8 +777,12 @@ static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
return;
}
+ /* Now, we can drift */
+
+ /* Get some information first */
const double timeBase = e->timeBase;
const int ti_old = c->ti_old;
+ const int ti_current = e->ti_current;
struct part *const parts = c->parts;
struct xpart *const xparts = c->xparts;
struct gpart *const gparts = c->gparts;
@@ -796,7 +795,7 @@ static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
if (!c->split) {
/* Check that we are actually going to move forward. */
- if (ti_current >= ti_old) {
+ if (ti_current > ti_old) {
/* Loop over all the g-particles in the cell */
const size_t nr_gparts = c->gcount;
@@ -840,6 +839,12 @@ static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
dx_max = sqrtf(dx2_max);
} /* Check that we are actually going to move forward. */
+
+ else {
+ /* ti_old == ti_current, just keep the current cell values. */
+ h_max = c->h_max;
+ dx_max = c->dx_max;
+ }
}
/* Otherwise, aggregate data from children. */
@@ -901,27 +906,22 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const double timeBase = e->timeBase;
- const int ti_current = r->e->ti_current;
const int count = c->count;
const int gcount = c->gcount;
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
struct gpart *restrict gparts = c->gparts;
- const double const_G = r->e->physical_constants->const_newton_G;
+ const double const_G = e->physical_constants->const_newton_G;
TIMER_TIC;
/* Anything to do here? */
- if (c->ti_end_min > ti_current) {
+ if (!cell_is_active(c, e)) {
c->updated = 0;
c->g_updated = 0;
return;
}
-#ifdef TASK_VERBOSE
- OUT;
-#endif
-
int updated = 0, g_updated = 0;
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
@@ -937,7 +937,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
/* If the g-particle has no counterpart and needs to be kicked */
if (gp->id_or_neg_offset > 0) {
- if (gp->ti_end <= ti_current) {
+ if (gpart_is_active(gp, e)) {
/* First, finish the force calculation */
gravity_end_force(gp, const_G);
@@ -968,7 +968,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xp = &xparts[k];
/* If particle needs to be kicked */
- if (p->ti_end <= ti_current) {
+ if (part_is_active(p, e)) {
/* First, finish the force loop */
hydro_end_force(p);
@@ -1126,7 +1126,7 @@ void *runner_main(void *data) {
/* Check that we haven't scheduled an inactive task */
#ifdef SWIFT_DEBUG_CHECKS
if (cj == NULL) { /* self */
- if (ci->ti_end_min > e->ti_current && t->type != task_type_sort)
+ if (!cell_is_active(ci, e) && t->type != task_type_sort)
error(
"Task (type='%s/%s') should have been skipped ti_current=%d "
"c->ti_end_min=%d",
@@ -1134,7 +1134,7 @@ void *runner_main(void *data) {
ci->ti_end_min);
/* Special case for sorts */
- if (ci->ti_end_min > e->ti_current && t->type == task_type_sort &&
+ if (!cell_is_active(ci, e) && t->type == task_type_sort &&
t->flags == 0)
error(
"Task (type='%s/%s') should have been skipped ti_current=%d "
@@ -1143,7 +1143,7 @@ void *runner_main(void *data) {
ci->ti_end_min, t->flags);
} else { /* pair */
- if (ci->ti_end_min > e->ti_current && cj->ti_end_min > e->ti_current)
+ if (!cell_is_active(ci, e) && !cell_is_active(cj, e))
error(
"Task (type='%s/%s') should have been skipped ti_current=%d "
"ci->ti_end_min=%d cj->ti_end_min=%d",
@@ -1224,11 +1224,11 @@ void *runner_main(void *data) {
runner_do_init(r, ci, 1);
break;
case task_type_ghost:
- runner_do_ghost(r, ci);
+ runner_do_ghost(r, ci, 1);
break;
#ifdef EXTRA_HYDRO_LOOP
case task_type_extra_ghost:
- runner_do_extra_ghost(r, ci);
+ runner_do_extra_ghost(r, ci, 1);
break;
#endif
case task_type_kick:
diff --git a/src/runner.h b/src/runner.h
index bffd3562df76d449d2bca762458c3d0d9b084b35..a8caf24248c99438f16729e2cac3e1031535f62b 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -48,11 +48,13 @@ struct runner {
};
/* Function prototypes. */
-void runner_do_ghost(struct runner *r, struct cell *c);
+void runner_do_ghost(struct runner *r, struct cell *c, int timer);
+void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer);
void runner_do_sort(struct runner *r, struct cell *c, int flag, int clock);
void runner_do_kick(struct runner *r, struct cell *c, int timer);
void runner_do_init(struct runner *r, struct cell *c, int timer);
void runner_do_cooling(struct runner *r, struct cell *c, int timer);
+void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
void *runner_main(void *data);
void runner_do_drift_mapper(void *map_data, int num_elements, void *extra_data);
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index 3c968cbf7d955198ad6bb44ab70e93af17735e99..c067d3bc9a576ee9c7dfb6e910eb1aa01012f755 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -109,7 +109,6 @@ void DOPAIR_NAIVE(struct runner *r, struct cell *restrict ci,
struct cell *restrict cj) {
const struct engine *e = r->e;
- const int ti_current = e->ti_current;
error("Don't use in actual runs ! Slow code !");
@@ -124,7 +123,7 @@ void DOPAIR_NAIVE(struct runner *r, struct cell *restrict ci,
TIMER_TIC;
/* Anything to do here? */
- if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
+ if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
const int count_i = ci->count;
const int count_j = cj->count;
@@ -210,7 +209,7 @@ void DOPAIR_NAIVE(struct runner *r, struct cell *restrict ci,
void DOSELF_NAIVE(struct runner *r, struct cell *restrict c) {
- const int ti_current = r->e->ti_current;
+ const struct engine *e = r->e;
error("Don't use in actual runs ! Slow code !");
@@ -226,7 +225,7 @@ void DOSELF_NAIVE(struct runner *r, struct cell *restrict c) {
TIMER_TIC;
/* Anything to do here? */
- if (c->ti_end_min > ti_current) return;
+ if (!cell_is_active(c, e)) return;
const int count = c->count;
struct part *restrict parts = c->parts;
@@ -706,8 +705,7 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
*/
void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
- struct engine *restrict e = r->e;
- const int ti_current = e->ti_current;
+ const struct engine *restrict e = r->e;
#ifdef WITH_VECTORIZATION
int icount = 0;
@@ -721,7 +719,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
TIMER_TIC;
/* Anything to do here? */
- if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
+ if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
/* Get the sort ID. */
double shift[3] = {0.0, 0.0, 0.0};
@@ -756,7 +754,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
/* Get a hold of the ith part in ci. */
struct part *restrict pi = &parts_i[sort_i[pid].i];
- if (pi->ti_end > ti_current) continue;
+ if (!part_is_active(pi, e)) continue;
const float hi = pi->h;
const double di = sort_i[pid].d + hi * kernel_gamma + dx_max - rshift;
if (di < dj_min) continue;
@@ -818,7 +816,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
/* Get a hold of the jth part in cj. */
struct part *restrict pj = &parts_j[sort_j[pjd].i];
- if (pj->ti_end > ti_current) continue;
+ if (!part_is_active(pj, e)) continue;
const float hj = pj->h;
const double dj = sort_j[pjd].d - hj * kernel_gamma - dx_max - rshift;
if (dj > di_max) continue;
@@ -894,7 +892,6 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
struct engine *restrict e = r->e;
- const int ti_current = e->ti_current;
#ifdef WITH_VECTORIZATION
int icount1 = 0;
@@ -914,7 +911,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
TIMER_TIC;
/* Anything to do here? */
- if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
+ if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
/* Get the shift ID. */
double shift[3] = {0.0, 0.0, 0.0};
@@ -946,28 +943,28 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
/* Collect the number of parts left and right below dt. */
int countdt_i = 0, countdt_j = 0;
struct entry *restrict sortdt_i = NULL, *restrict sortdt_j = NULL;
- if (ci->ti_end_max <= ti_current) {
+ if (cell_is_all_active(ci, e)) {
sortdt_i = sort_i;
countdt_i = count_i;
- } else if (ci->ti_end_min <= ti_current) {
+ } else if (cell_is_active(ci, e)) {
if (posix_memalign((void *)&sortdt_i, VEC_SIZE * sizeof(float),
sizeof(struct entry) * count_i) != 0)
error("Failed to allocate dt sortlists.");
for (int k = 0; k < count_i; k++)
- if (parts_i[sort_i[k].i].ti_end <= ti_current) {
+ if (part_is_active(&parts_i[sort_i[k].i], e)) {
sortdt_i[countdt_i] = sort_i[k];
countdt_i += 1;
}
}
- if (cj->ti_end_max <= ti_current) {
+ if (cell_is_all_active(cj, e)) {
sortdt_j = sort_j;
countdt_j = count_j;
- } else if (cj->ti_end_min <= ti_current) {
+ } else if (cell_is_active(cj, e)) {
if (posix_memalign((void *)&sortdt_j, VEC_SIZE * sizeof(float),
sizeof(struct entry) * count_j) != 0)
error("Failed to allocate dt sortlists.");
for (int k = 0; k < count_j; k++)
- if (parts_j[sort_j[k].i].ti_end <= ti_current) {
+ if (part_is_active(&parts_j[sort_j[k].i], e)) {
sortdt_j[countdt_j] = sort_j[k];
countdt_j += 1;
}
@@ -988,7 +985,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
const float hig2 = hi * hi * kernel_gamma2;
/* Look at valid dt parts only? */
- if (pi->ti_end > ti_current) {
+ if (!part_is_active(pi, e)) {
/* Loop over the parts in cj within dt. */
for (int pjd = 0; pjd < countdt_j && sortdt_j[pjd].d < di; pjd++) {
@@ -1062,7 +1059,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
#ifndef WITH_VECTORIZATION
/* Does pj need to be updated too? */
- if (pj->ti_end <= ti_current)
+ if (part_is_active(pj, e))
IACT(r2, dx, hi, hj, pi, pj);
else
IACT_NONSYM(r2, dx, hi, hj, pi, pj);
@@ -1070,7 +1067,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
#else
/* Does pj need to be updated too? */
- if (pj->ti_end <= ti_current) {
+ if (part_is_active(pj, e)) {
/* Add this interaction to the symmetric queue. */
r2q2[icount2] = r2;
@@ -1132,7 +1129,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
const float hjg2 = hj * hj * kernel_gamma2;
/* Is this particle outside the dt? */
- if (pj->ti_end > ti_current) {
+ if (!part_is_active(pj, e)) {
/* Loop over the parts in ci. */
for (int pid = countdt_i - 1; pid >= 0 && sortdt_i[pid].d > dj; pid--) {
@@ -1205,7 +1202,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
#ifndef WITH_VECTORIZATION
/* Does pi need to be updated too? */
- if (pi->ti_end <= ti_current)
+ if (part_is_active(pi, e))
IACT(r2, dx, hj, hi, pj, pi);
else
IACT_NONSYM(r2, dx, hj, hi, pj, pi);
@@ -1213,7 +1210,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
#else
/* Does pi need to be updated too? */
- if (pi->ti_end <= ti_current) {
+ if (part_is_active(pi, e)) {
/* Add this interaction to the symmetric queue. */
r2q2[icount2] = r2;
@@ -1270,10 +1267,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
IACT(r2q2[k], &dxq2[3 * k], hiq2[k], hjq2[k], piq2[k], pjq2[k]);
#endif
- if (ci->ti_end_max > ti_current && ci->ti_end_min <= ti_current)
- free(sortdt_i);
- if (cj->ti_end_max > ti_current && cj->ti_end_min <= ti_current)
- free(sortdt_j);
+ /* Clean-up if necessary */
+ if (cell_is_active(ci, e) && !cell_is_all_active(ci, e)) free(sortdt_i);
+ if (cell_is_active(cj, e) && !cell_is_all_active(cj, e)) free(sortdt_j);
TIMER_TOC(TIMER_DOPAIR);
}
@@ -1286,7 +1282,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
*/
void DOSELF1(struct runner *r, struct cell *restrict c) {
- const int ti_current = r->e->ti_current;
+ const struct engine *e = r->e;
#ifdef WITH_VECTORIZATION
int icount1 = 0;
@@ -1305,8 +1301,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
TIMER_TIC;
- if (c->ti_end_min > ti_current) return;
- if (c->ti_end_max < ti_current) error("Cell in an impossible time-zone");
+ if (!cell_is_active(c, e)) return;
struct part *restrict parts = c->parts;
const int count = c->count;
@@ -1318,7 +1313,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
count * sizeof(int)) != 0)
error("Failed to allocate indt.");
for (int k = 0; k < count; k++)
- if (parts[k].ti_end <= ti_current) {
+ if (part_is_active(&parts[k], e)) {
indt[countdt] = k;
countdt += 1;
}
@@ -1336,7 +1331,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
const float hig2 = hi * hi * kernel_gamma2;
/* Is the ith particle inactive? */
- if (pi->ti_end > ti_current) {
+ if (!part_is_active(pi, e)) {
/* Loop over the other particles .*/
for (int pjd = firstdt; pjd < countdt; pjd++) {
@@ -1407,7 +1402,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
r2 += dx[k] * dx[k];
}
const int doj =
- (pj->ti_end <= ti_current) && (r2 < hj * hj * kernel_gamma2);
+ (part_is_active(pj, e)) && (r2 < hj * hj * kernel_gamma2);
/* Hit or miss? */
if (r2 < hig2 || doj) {
@@ -1518,7 +1513,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
*/
void DOSELF2(struct runner *r, struct cell *restrict c) {
- const int ti_current = r->e->ti_current;
+ const struct engine *e = r->e;
#ifdef WITH_VECTORIZATION
int icount1 = 0;
@@ -1537,8 +1532,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
TIMER_TIC;
- if (c->ti_end_min > ti_current) return;
- if (c->ti_end_max < ti_current) error("Cell in an impossible time-zone");
+ if (!cell_is_active(c, e)) return;
struct part *restrict parts = c->parts;
const int count = c->count;
@@ -1550,7 +1544,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
count * sizeof(int)) != 0)
error("Failed to allocate indt.");
for (int k = 0; k < count; k++)
- if (parts[k].ti_end <= ti_current) {
+ if (part_is_active(&parts[k], e)) {
indt[countdt] = k;
countdt += 1;
}
@@ -1568,7 +1562,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
const float hig2 = hi * hi * kernel_gamma2;
/* Is the ith particle not active? */
- if (pi->ti_end > ti_current) {
+ if (!part_is_active(pi, e)) {
/* Loop over the other particles .*/
for (int pjd = firstdt; pjd < countdt; pjd++) {
@@ -1645,7 +1639,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
#ifndef WITH_VECTORIZATION
/* Does pj need to be updated too? */
- if (pj->ti_end <= ti_current)
+ if (part_is_active(pj, e))
IACT(r2, dx, hi, hj, pi, pj);
else
IACT_NONSYM(r2, dx, hi, hj, pi, pj);
@@ -1653,7 +1647,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
#else
/* Does pj need to be updated too? */
- if (pj->ti_end <= ti_current) {
+ if (part_is_active(pj, e)) {
/* Add this interaction to the symmetric queue. */
r2q2[icount2] = r2;
@@ -1731,12 +1725,12 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid,
int gettimer) {
struct space *s = r->e->s;
- const int ti_current = r->e->ti_current;
+ const struct engine *e = r->e;
TIMER_TIC;
/* Should we even bother? */
- if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
+ if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
/* Get the cell dimensions. */
const float h = min(ci->width[0], min(ci->width[1], ci->width[2]));
@@ -1950,7 +1944,7 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid,
}
/* Otherwise, compute the pair directly. */
- else if (ci->ti_end_min <= ti_current || cj->ti_end_min <= ti_current) {
+ else if (cell_is_active(ci, e) || cell_is_active(cj, e)) {
/* Do any of the cells need to be sorted first? */
if (!(ci->sorted & (1 << sid))) runner_do_sort(r, ci, (1 << sid), 1);
@@ -1972,12 +1966,10 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid,
*/
void DOSUB_SELF1(struct runner *r, struct cell *ci, int gettimer) {
- const int ti_current = r->e->ti_current;
-
TIMER_TIC;
/* Should we even bother? */
- if (ci->ti_end_min > ti_current) return;
+ if (!cell_is_active(ci, r->e)) return;
/* Recurse? */
if (ci->split) {
@@ -2014,13 +2006,13 @@ void DOSUB_SELF1(struct runner *r, struct cell *ci, int gettimer) {
void DOSUB_PAIR2(struct runner *r, struct cell *ci, struct cell *cj, int sid,
int gettimer) {
- struct space *s = r->e->s;
- const int ti_current = r->e->ti_current;
+ const struct engine *e = r->e;
+ struct space *s = e->s;
TIMER_TIC;
/* Should we even bother? */
- if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
+ if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
/* Get the cell dimensions. */
const float h = min(ci->width[0], min(ci->width[1], ci->width[2]));
@@ -2234,7 +2226,7 @@ void DOSUB_PAIR2(struct runner *r, struct cell *ci, struct cell *cj, int sid,
}
/* Otherwise, compute the pair directly. */
- else if (ci->ti_end_min <= ti_current || cj->ti_end_min <= ti_current) {
+ else if (cell_is_active(ci, e) || cell_is_active(cj, e)) {
/* Do any of the cells need to be sorted first? */
if (!(ci->sorted & (1 << sid))) runner_do_sort(r, ci, (1 << sid), 1);
@@ -2256,12 +2248,10 @@ void DOSUB_PAIR2(struct runner *r, struct cell *ci, struct cell *cj, int sid,
*/
void DOSUB_SELF2(struct runner *r, struct cell *ci, int gettimer) {
- const int ti_current = r->e->ti_current;
-
TIMER_TIC;
/* Should we even bother? */
- if (ci->ti_end_min > ti_current) return;
+ if (!cell_is_active(ci, r->e)) return;
/* Recurse? */
if (ci->split) {
@@ -2287,8 +2277,8 @@ void DOSUB_SELF2(struct runner *r, struct cell *ci, int gettimer) {
void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
int *ind, int count, struct cell *cj, int sid, int gettimer) {
- struct space *s = r->e->s;
- const int ti_current = r->e->ti_current;
+ const struct engine *e = r->e;
+ struct space *s = e->s;
TIMER_TIC;
@@ -2850,7 +2840,7 @@ void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
}
/* Otherwise, compute the pair directly. */
- else if (ci->ti_end_min <= ti_current || cj->ti_end_min <= ti_current) {
+ else if (cell_is_active(ci, e) || cell_is_active(cj, e)) {
/* Get the relative distance between the pairs, wrapping. */
double shift[3] = {0.0, 0.0, 0.0};
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 0fcd2d2e80a72b92588acd5b8275b9dafc68df45..59a5ae496680390c23458bde65b4bba321ffe7a1 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -71,7 +71,6 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
const int gcount = ci->gcount;
struct gpart *restrict gparts = ci->gparts;
const struct multipole multi = cj->multipole;
- const int ti_current = e->ti_current;
const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
TIMER_TIC;
@@ -84,7 +83,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
#endif
/* Anything to do here? */
- if (ci->ti_end_min > ti_current) return;
+ if (!cell_is_active(ci, e)) return;
#if ICHECK > 0
for (int pid = 0; pid < gcount; pid++) {
@@ -105,7 +104,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
/* Get a hold of the ith part in ci. */
struct gpart *restrict gp = &gparts[pid];
- if (gp->ti_end > ti_current) continue;
+ if (!gpart_is_active(gp, e)) continue;
/* Compute the pairwise distance. */
const float dx[3] = {multi.CoM[0] - gp->x[0], // x
@@ -138,7 +137,6 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
const int gcount_j = cj->gcount;
struct gpart *restrict gparts_i = ci->gparts;
struct gpart *restrict gparts_j = cj->gparts;
- const int ti_current = e->ti_current;
const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
TIMER_TIC;
@@ -150,7 +148,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
#endif
/* Anything to do here? */
- if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
+ if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
#if ICHECK > 0
for (int pid = 0; pid < gcount_i; pid++) {
@@ -195,17 +193,17 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Interact ! */
- if (gpi->ti_end <= ti_current && gpj->ti_end <= ti_current) {
+ if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) {
runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
} else {
- if (gpi->ti_end <= ti_current) {
+ if (gpart_is_active(gpi, e)) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
- } else if (gpj->ti_end <= ti_current) {
+ } else if (gpart_is_active(gpj, e)) {
dx[0] = -dx[0];
dx[1] = -dx[1];
@@ -233,7 +231,6 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
const struct engine *e = r->e;
const int gcount = c->gcount;
struct gpart *restrict gparts = c->gparts;
- const int ti_current = e->ti_current;
const float rlr_inv = 1. / (const_gravity_a_smooth * c->super->width[0]);
TIMER_TIC;
@@ -244,7 +241,7 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
#endif
/* Anything to do here? */
- if (c->ti_end_min > ti_current) return;
+ if (!cell_is_active(c, e)) return;
#if ICHECK > 0
for (int pid = 0; pid < gcount; pid++) {
@@ -278,17 +275,17 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Interact ! */
- if (gpi->ti_end <= ti_current && gpj->ti_end <= ti_current) {
+ if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) {
runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
} else {
- if (gpi->ti_end <= ti_current) {
+ if (gpart_is_active(gpi, e)) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
- } else if (gpj->ti_end <= ti_current) {
+ } else if (gpart_is_active(gpj, e)) {
dx[0] = -dx[0];
dx[1] = -dx[1];
@@ -490,14 +487,13 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
const struct engine *e = r->e;
struct cell *cells = e->s->cells_top;
const int nr_cells = e->s->nr_cells;
- const int ti_current = e->ti_current;
const double max_d =
const_gravity_a_smooth * const_gravity_r_cut * ci->width[0];
const double max_d2 = max_d * max_d;
const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]};
/* Anything to do here? */
- if (ci->ti_end_min > ti_current) return;
+ if (!cell_is_active(ci, e)) return;
/* Loop over all the cells and go for a p-m interaction if far enough but not
* too far */
diff --git a/src/scheduler.c b/src/scheduler.c
index 593db887434301eef8f43460b24d403d71954bba..16e513a4189aca6615e545d53c1773d07209329f 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1387,10 +1387,11 @@ struct task *scheduler_gettask(struct scheduler *s, int qid,
/* If we failed, take a short nap. */
#ifdef WITH_MPI
- if (res == NULL && qid > 1) {
+ if (res == NULL && qid > 1)
#else
- if (res == NULL) {
+ if (res == NULL)
#endif
+ {
pthread_mutex_lock(&s->sleep_mutex);
res = queue_gettask(&s->queues[qid], prev, 1);
if (res == NULL && s->waiting > 0) {
diff --git a/src/space.c b/src/space.c
index 18f0e21e87b22715e6318be5dfbc3c0cc9a2cc12..5f7e3522eef8c398faba1aee4b03c49ab94232b1 100644
--- a/src/space.c
+++ b/src/space.c
@@ -346,7 +346,6 @@ void space_regrid(struct space *s, int verbose) {
if (verbose)
message("set cell dimensions to [ %i %i %i ].", cdim[0], cdim[1],
cdim[2]);
- fflush(stdout);
#ifdef WITH_MPI
if (oldnodeIDs != NULL) {
@@ -1563,6 +1562,15 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
else
c->owner = 0; /* Ok, there is really nothing on this rank... */
}
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* All cells and particles should have consistent h_max values. */
+ for (int ind = 0; ind < num_cells; ind++) {
+ int depth = 0;
+ if (!checkCellhdxmax(&cells_top[ind], &depth))
+ message(" at cell depth %d", depth);
+ }
+#endif
}
/**
diff --git a/src/timers.h b/src/timers.h
index d75508afd88cf2fd57d06f9530bff8607d79d127..bc877d4094425a4948290d2c7c099f49cbd44280 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -56,6 +56,7 @@ enum {
timer_dosub_pair_grav,
timer_dopair_subset,
timer_do_ghost,
+ timer_do_extra_ghost,
timer_dorecv_cell,
timer_gettask,
timer_qget,
diff --git a/src/tools.c b/src/tools.c
index 060bf1439f30dc6237938c060bc4ddc8d9be822b..e526bb1b838f6d97b72eadb4070f3f2a94938c04 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -481,7 +481,7 @@ void engine_single_density(double *dim, long long int pid,
}
/* Dump the result. */
- hydro_end_density(&p, 0);
+ hydro_end_density(&p);
message("part %lli (h=%e) has wcount=%e, rho=%e.", p.id, p.h,
p.density.wcount, hydro_get_density(&p));
fflush(stdout);
diff --git a/tests/test125cells.c b/tests/test125cells.c
index d423efba3fda764310721586e14192afdee18a85..bb92d8e4b0fb8eb19afefc268a03ac0faa93737b 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -617,7 +617,7 @@ int main(int argc, char *argv[]) {
#endif
/* Ghost to finish everything on the central cells */
- for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j]);
+ for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j], 0);
/* Do the force calculation */
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
@@ -704,7 +704,7 @@ int main(int argc, char *argv[]) {
#endif
/* Ghost to finish everything on the central cells */
- for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j]);
+ for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j], 0);
/* Do the force calculation */
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 22619af53c39218da34d771fab6ed2d10993689c..f58b4dc410637f3d91369dab1b442de0b7044c08 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -167,7 +167,7 @@ void zero_particle_fields(struct cell *c) {
*/
void end_calculation(struct cell *c) {
for (int pid = 0; pid < c->count; pid++) {
- hydro_end_density(&c->parts[pid], 1);
+ hydro_end_density(&c->parts[pid]);
}
}