diff --git a/.gitignore b/.gitignore
index f85ebc57eb97c2e778fd57ad701d9b907168c526..798b187d4b326fb24bd4ccccabf7dc84b9bb20fd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -30,6 +30,8 @@ examples/*/*/*.mp4
examples/*/*/*.txt
examples/*/*/*.rst
examples/*/*/*.hdf5
+examples/*/*/*.csv
+examples/*/*/*.dot
examples/*/*/restart/*
examples/*/*/used_parameters.yml
examples/*/*/log*
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index fe48398e4b0dd4bfdede6781e708bd222b6eedae..db841f347f681c42c1b305c2d130ee0b55d639ae 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -775,6 +775,7 @@ INPUT += @top_srcdir@/src/chemistry/EAGLE
INPUT += @top_srcdir@/src/entropy_floor/EAGLE
INPUT += @top_srcdir@/src/star_formation/EAGLE
INPUT += @top_srcdir@/src/tracers/EAGLE
+INPUT += @top_srcdir@/src/stars/EAGLE
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index 90070f6a342d60d48d8fd2c794014f9d638030be..d09bbb51e90d843dd6731c5fcbd48b9c586713f9 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -57,7 +57,7 @@ InitialConditions:
periodic: 1
cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
cleanup_velocity_factors: 1 # Remove the sqrt(a) factor in the velocities inherited from Gadget
-
+
EAGLEChemistry: # Solar abundances
init_abundance_metal: 0.014
init_abundance_Hydrogen: 0.70649785
@@ -76,3 +76,31 @@ EAGLECooling:
He_reion_z_centre: 3.5
He_reion_z_sigma: 0.5
He_reion_eV_p_H: 2.0
+
+# EAGLE star formation parameters
+EAGLEStarFormation:
+ EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
+ EOS_temperature_norm_K: 8000 # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
+ EOS_gamma_effective: 1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
+ gas_fraction: 0.3 # The gas fraction used internally by the model.
+ KS_normalisation: 1.515e-4 # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+ KS_exponent: 1.4 # The exponent of the Kennicutt-Schmidt law.
+ KS_min_over_density: 57.7 # The over-density above which star-formation is allowed.
+ KS_high_density_threshold_H_p_cm3: 1e3 # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+ KS_high_density_exponent: 2.0 # Slope of the Kennicut-Schmidt law above the high-density threshold.
+ KS_temperature_margin_dex: 0.5 # Logarithm base 10 of the maximal temperature difference above the EOS allowed to form stars.
+ threshold_norm_H_p_cm3: 0.1 # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+ threshold_Z0: 0.002 # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+ threshold_slope: -0.64 # Slope of the metal-dependant star formation threshold
+ threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+
+# Parameters for the EAGLE "equation of state"
+EAGLEEntropyFloor:
+ Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+ Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+ Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+ Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+ Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+ Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index 3d7d5de4fe2c98a04c986840bc2804edaf780c17..75799647b4e95ebd75202748c67a2c18c423f532 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -84,3 +84,30 @@ EAGLECooling:
He_reion_z_sigma: 0.5
He_reion_eV_p_H: 2.0
+# EAGLE star formation parameters
+EAGLEStarFormation:
+ EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
+ EOS_temperature_norm_K: 8000 # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
+ EOS_gamma_effective: 1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
+ gas_fraction: 0.3 # The gas fraction used internally by the model.
+ KS_normalisation: 1.515e-4 # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+ KS_exponent: 1.4 # The exponent of the Kennicutt-Schmidt law.
+ KS_min_over_density: 57.7 # The over-density above which star-formation is allowed.
+ KS_high_density_threshold_H_p_cm3: 1e3 # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+ KS_high_density_exponent: 2.0 # Slope of the Kennicut-Schmidt law above the high-density threshold.
+ KS_temperature_margin_dex: 0.5 # Logarithm base 10 of the maximal temperature difference above the EOS allowed to form stars.
+ threshold_norm_H_p_cm3: 0.1 # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+ threshold_Z0: 0.002 # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+ threshold_slope: -0.64 # Slope of the metal-dependant star formation threshold
+ threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+
+# Parameters for the EAGLE "equation of state"
+EAGLEEntropyFloor:
+ Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+ Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+ Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+ Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+ Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+ Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index 8d622998c5ea1c754e877643e42ee3d97bfa720e..6c0c7421ba4f804437a8086b42fb2878bd3904b1 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -77,4 +77,32 @@ EAGLECooling:
H_reion_z: 11.5
He_reion_z_centre: 3.5
He_reion_z_sigma: 0.5
- He_reion_eV_p_H: 2.0
\ No newline at end of file
+ He_reion_eV_p_H: 2.0
+
+# EAGLE star formation parameters
+EAGLEStarFormation:
+ EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
+ EOS_temperature_norm_K: 8000 # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
+ EOS_gamma_effective: 1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
+ gas_fraction: 0.3 # The gas fraction used internally by the model.
+ KS_normalisation: 1.515e-4 # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+ KS_exponent: 1.4 # The exponent of the Kennicutt-Schmidt law.
+ KS_min_over_density: 57.7 # The over-density above which star-formation is allowed.
+ KS_high_density_threshold_H_p_cm3: 1e3 # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+ KS_high_density_exponent: 2.0 # Slope of the Kennicut-Schmidt law above the high-density threshold.
+ KS_temperature_margin_dex: 0.5 # Logarithm base 10 of the maximal temperature difference above the EOS allowed to form stars.
+ threshold_norm_H_p_cm3: 0.1 # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+ threshold_Z0: 0.002 # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+ threshold_slope: -0.64 # Slope of the metal-dependant star formation threshold
+ threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+
+# Parameters for the EAGLE "equation of state"
+EAGLEEntropyFloor:
+ Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+ Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+ Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+ Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+ Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+ Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index fb885fa44f119bb3e296947997f71581d5ca0f48..313a5a384324eecd56455ea22bbc96d147982d6b 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -87,3 +87,30 @@ EAGLECooling:
He_reion_z_centre: 3.5
He_reion_z_sigma: 0.5
He_reion_eV_p_H: 2.0
+# EAGLE star formation parameters
+EAGLEStarFormation:
+ EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
+ EOS_temperature_norm_K: 8000 # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
+ EOS_gamma_effective: 1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
+ gas_fraction: 0.3 # The gas fraction used internally by the model.
+ KS_normalisation: 1.515e-4 # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+ KS_exponent: 1.4 # The exponent of the Kennicutt-Schmidt law.
+ KS_min_over_density: 57.7 # The over-density above which star-formation is allowed.
+ KS_high_density_threshold_H_p_cm3: 1e3 # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+ KS_high_density_exponent: 2.0 # Slope of the Kennicut-Schmidt law above the high-density threshold.
+ KS_temperature_margin_dex: 0.5 # Logarithm base 10 of the maximal temperature difference above the EOS allowed to form stars.
+ threshold_norm_H_p_cm3: 0.1 # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+ threshold_Z0: 0.002 # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+ threshold_slope: -0.64 # Slope of the metal-dependant star formation threshold
+ threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+
+# Parameters for the EAGLE "equation of state"
+EAGLEEntropyFloor:
+ Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+ Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+ Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+ Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+ Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+ Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
diff --git a/examples/GEAR/DwarfGalaxy/dwarf_galaxy.yml b/examples/GEAR/DwarfGalaxy/dwarf_galaxy.yml
index 4c5e2a82b017725929138de011b1f3ed1fe9f1ef..00fd889d4f58a4dadccc52ef5c6d8315ac2a4012 100644
--- a/examples/GEAR/DwarfGalaxy/dwarf_galaxy.yml
+++ b/examples/GEAR/DwarfGalaxy/dwarf_galaxy.yml
@@ -32,7 +32,7 @@ TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 1. # The end time of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
+ dt_max: 1e-5 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml
index 2a885f5770f88c3ee6dd75e194d0fa713160e39e..7ba1e601c764d9c12b93178efd8226601af8373c 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml
@@ -36,8 +36,9 @@ Statistics:
# Parameters related to the initial conditions
InitialConditions:
file_name: fid.hdf5 # The file to read
- periodic: 0 # Are we running with periodic ICs?
-
+ periodic: 0 # Are we running with periodic ICs?
+ stars_smoothing_length: 0.5
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py
index 4e5fcf2ceb2b3f0682871cc533fab93039ff38e9..89a87923148cb6872ab17b6d7229aef597ef3287 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py
@@ -101,7 +101,7 @@ gas_sSFR = gas_SFR / gas_mass
# Read the Star properties
stars_pos = f["/PartType4/Coordinates"][:, :]
stars_BirthDensity = f["/PartType4/BirthDensity"][:]
-stars_BirthTime = f["/PartType4/Birth_time"][:]
+stars_BirthTime = f["/PartType4/BirthTime"][:]
stars_XH = f["/PartType4/ElementAbundance"][:,0]
# Centre the box
diff --git a/examples/main.c b/examples/main.c
index 1d46ee954b55c5965a6ab65dae5b5038e4c72394..6fc5b433719822558d531f4ed2691e7127139a79 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -347,7 +347,17 @@ int main(int argc, char *argv[]) {
if (myrank == 0) {
argparse_usage(&argparse);
printf(
- "\nError: Cannot process feedback without stars, -S must be "
+ "\nError: Cannot process feedback without stars, --stars must be "
+ "chosen.\n");
+ }
+ return 1;
+ }
+
+ if (!with_hydro && with_feedback) {
+ if (myrank == 0) {
+ argparse_usage(&argparse);
+ printf(
+ "\nError: Cannot process feedback without gas, --hydro must be "
"chosen.\n");
}
return 1;
@@ -460,7 +470,6 @@ int main(int argc, char *argv[]) {
#ifdef WITH_MPI
if (with_mpole_reconstruction && nr_nodes > 1)
error("Cannot reconstruct m-poles every step over MPI (yet).");
- if (with_feedback) error("Can't run with feedback over MPI (yet).");
if (with_star_formation)
error("Can't run with star formation over MPI (yet)");
if (with_limiter) error("Can't run with time-step limiter over MPI (yet)");
@@ -893,10 +902,11 @@ int main(int argc, char *argv[]) {
if (myrank == 0) cooling_print(&cooling_func);
/* Initialise the star formation law and its properties */
+ bzero(&starform, sizeof(struct star_formation));
if (with_star_formation)
starformation_init(params, &prog_const, &us, &hydro_properties,
&starform);
- if (myrank == 0) starformation_print(&starform);
+ if (with_star_formation && myrank == 0) starformation_print(&starform);
/* Initialise the chemistry */
bzero(&chemistry, sizeof(struct chemistry_global_data));
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 18a6a4778c98ed73e1674908a9bbedbf3e38f013..22bbf3db4f4f49f1cce6c1aa817b8228f829437f 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -65,8 +65,6 @@ Scheduler:
cell_sub_size_self_hydro: 32000 # (Optional) Maximal number of interactions per sub-self hydro task (this is the default value).
cell_sub_size_pair_grav: 256000000 # (Optional) Maximal number of interactions per sub-pair gravity task (this is the default value).
cell_sub_size_self_grav: 32000 # (Optional) Maximal number of interactions per sub-self gravity task (this is the default value).
- cell_sub_size_pair_stars: 256000000 # (Optional) Maximal number of interactions per sub-pair stars task (this is the default value).
- cell_sub_size_self_stars: 32000 # (Optional) Maximal number of interactions per sub-self stars task (this is the default value).
cell_split_size: 400 # (Optional) Maximal number of particles per cell (this is the default value).
cell_subdepth_diff_grav: 4 # (Optional) Maximal depth difference between leaves and a cell that gravity tasks can be pushed down to (this is the default value).
cell_extra_parts: 0 # (Optional) Number of spare parts per top-level allocated at rebuild time for on-the-fly creation.
diff --git a/src/active.h b/src/active.h
index 5bbbd3803cb09e7aa05ddb15e2e5c2a15b27602c..6466cd314fdc18ad324bf01a1ff4e73e214e35d5 100644
--- a/src/active.h
+++ b/src/active.h
@@ -85,9 +85,16 @@ __attribute__((always_inline)) INLINE static int cell_are_gpart_drifted(
__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);
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->stars.ti_old_part > e->ti_current)
+ error(
+ "Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
+ "and e->ti_current=%lld (t=%e)",
+ c->stars.ti_old_part, c->stars.ti_old_part * e->time_base,
+ e->ti_current, e->ti_current * e->time_base);
+#endif
+
+ return (c->stars.ti_old_part == e->ti_current);
}
/* Are cells / particles active for regular tasks ? */
@@ -125,7 +132,7 @@ __attribute__((always_inline)) INLINE static int cell_is_all_active_hydro(
const struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
- if (c->hydro.ti_end_max < e->ti_current)
+ if (c->hydro.count > 0 && c->hydro.ti_end_max < e->ti_current)
error(
"cell in an impossible time-zone! c->ti_end_max=%lld "
"e->ti_current=%lld",
@@ -181,7 +188,7 @@ __attribute__((always_inline)) INLINE static int cell_is_all_active_gravity(
const struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
- if (c->grav.ti_end_max < e->ti_current)
+ if (c->grav.count > 0 && c->grav.ti_end_max < e->ti_current)
error(
"cell in an impossible time-zone! c->ti_end_max=%lld "
"e->ti_current=%lld",
@@ -383,6 +390,28 @@ __attribute__((always_inline)) INLINE static int cell_is_starting_gravity(
return (c->grav.ti_beg_max == e->ti_current);
}
+/**
+ * @brief Does a cell contain any s-particle starting their time-step now ?
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #cell contains at least an active particle, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int cell_is_starting_stars(
+ const struct cell *c, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->stars.ti_beg_max > e->ti_current)
+ error(
+ "cell in an impossible time-zone! c->ti_beg_max=%lld (t=%e) and "
+ "e->ti_current=%lld (t=%e, a=%e)",
+ c->stars.ti_beg_max, c->stars.ti_beg_max * e->time_base, e->ti_current,
+ e->ti_current * e->time_base, e->cosmology->a);
+#endif
+
+ return (c->stars.ti_beg_max == e->ti_current);
+}
+
/**
* @brief Is this particle starting its time-step now ?
*
diff --git a/src/cache.h b/src/cache.h
index b6735f381da5c6e615d50e0c5768ca201022aa59..e5a62f33b3eb492f9da6e0e98ed767d6b8de32dd 100644
--- a/src/cache.h
+++ b/src/cache.h
@@ -123,8 +123,7 @@ __attribute__((always_inline)) INLINE void cache_init(struct cache *c,
size_t count) {
/* Align cache on correct byte boundary and pad cache size to be a multiple of
- * the vector size
- * and include 2 vector lengths for remainder operations. */
+ * the vector size and include 2 vector lengths for remainder operations. */
size_t pad = 2 * VEC_SIZE, rem = count % VEC_SIZE;
if (rem > 0) pad += VEC_SIZE - rem;
size_t sizeBytes = (count + pad) * sizeof(float);
diff --git a/src/cell.c b/src/cell.c
index 4f44abd924be97e54e52fffda8559f7a0dcfe68a..af4edc6a92dd8ed84dd4cd8cc869b8c32c4f81e8 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -381,14 +381,17 @@ int cell_pack(struct cell *restrict c, struct pcell *restrict pc,
/* Start by packing the data of the current cell. */
pc->hydro.h_max = c->hydro.h_max;
+ pc->stars.h_max = c->stars.h_max;
pc->hydro.ti_end_min = c->hydro.ti_end_min;
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->stars.ti_end_max = c->stars.ti_end_max;
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;
+ pc->stars.ti_old_part = c->stars.ti_old_part;
pc->hydro.count = c->hydro.count;
pc->grav.count = c->grav.count;
pc->stars.count = c->stars.count;
@@ -485,14 +488,17 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
/* Unpack the current pcell. */
c->hydro.h_max = pc->hydro.h_max;
+ c->stars.h_max = pc->stars.h_max;
c->hydro.ti_end_min = pc->hydro.ti_end_min;
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->stars.ti_end_max = pc->stars.ti_end_max;
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;
+ c->stars.ti_old_part = pc->stars.ti_old_part;
c->hydro.count = pc->hydro.count;
c->grav.count = pc->grav.count;
c->stars.count = pc->stars.count;
@@ -619,6 +625,7 @@ int cell_pack_end_step(struct cell *restrict c,
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].stars.ti_end_max = c->stars.ti_end_max;
pcells[0].hydro.dx_max_part = c->hydro.dx_max_part;
pcells[0].stars.dx_max_part = c->stars.dx_max_part;
@@ -657,6 +664,7 @@ int cell_unpack_end_step(struct cell *restrict c,
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->stars.ti_end_max = pcells[0].stars.ti_end_max;
c->hydro.dx_max_part = pcells[0].hydro.dx_max_part;
c->stars.dx_max_part = pcells[0].stars.dx_max_part;
@@ -1475,7 +1483,7 @@ void cell_check_part_drift_point(struct cell *c, void *data) {
}
/**
- * @brief Checks that the #gpart and #spart in a cell are at the
+ * @brief Checks that the #gpart in a cell are at the
* current point in time
*
* Calls error() if the cell is not at the current time.
@@ -1506,11 +1514,42 @@ void cell_check_gpart_drift_point(struct cell *c, void *data) {
c->grav.parts[i].time_bin != time_bin_inhibited)
error("g-part in an incorrect time-zone! gp->ti_drift=%lld ti_drift=%lld",
c->grav.parts[i].ti_drift, ti_drift);
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
+/**
+ * @brief Checks that the #spart in a cell are at the
+ * current point in time
+ *
+ * Calls error() if the cell is not at the current time.
+ *
+ * @param c Cell to act upon
+ * @param data The current time on the integer time-line
+ */
+void cell_check_spart_drift_point(struct cell *c, void *data) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ const integertime_t ti_drift = *(integertime_t *)data;
+
+ /* Only check local cells */
+ if (c->nodeID != engine_rank) return;
+
+ /* Only check cells with content */
+ if (c->stars.count == 0) return;
+
+ if (c->stars.ti_old_part != ti_drift)
+ error(
+ "Cell in an incorrect time-zone! c->stars.ti_old_part=%lld "
+ "ti_drift=%lld",
+ c->stars.ti_old_part, ti_drift);
for (int i = 0; i < c->stars.count; ++i)
if (c->stars.parts[i].ti_drift != ti_drift &&
c->stars.parts[i].time_bin != time_bin_inhibited)
- error("s-part in an incorrect time-zone! sp->ti_drift=%lld ti_drift=%lld",
+ error("g-part in an incorrect time-zone! gp->ti_drift=%lld ti_drift=%lld",
c->stars.parts[i].ti_drift, ti_drift);
#else
error("Calling debugging code without debugging flag activated.");
@@ -1800,6 +1839,8 @@ void cell_clear_drift_flags(struct cell *c, void *data) {
c->hydro.do_sub_drift = 0;
c->grav.do_drift = 0;
c->grav.do_sub_drift = 0;
+ c->stars.do_drift = 0;
+ c->stars.do_sub_drift = 0;
}
/**
@@ -1896,8 +1937,39 @@ 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);
+
+ /* If this cell is already marked for drift, quit early. */
+ if (c->stars.do_drift) return;
+
+ /* Mark this cell for drifting. */
+ c->stars.do_drift = 1;
+
+ /* Set the do_stars_sub_drifts all the way up and activate the super drift
+ if this has not yet been done. */
+ if (c == c->hydro.super) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->stars.drift == NULL)
+ error("Trying to activate un-existing c->stars.drift");
+#endif
+ scheduler_activate(s, c->stars.drift);
+ } else {
+ for (struct cell *parent = c->parent;
+ parent != NULL && !parent->stars.do_sub_drift;
+ parent = parent->parent) {
+
+ /* Mark this cell for drifting */
+ parent->stars.do_sub_drift = 1;
+
+ if (parent == c->hydro.super) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (parent->stars.drift == NULL)
+ error("Trying to activate un-existing parent->stars.drift");
+#endif
+ scheduler_activate(s, parent->stars.drift);
+ break;
+ }
+ }
+ }
}
/**
@@ -1905,7 +1977,7 @@ void cell_activate_drift_spart(struct cell *c, struct scheduler *s) {
*/
void cell_activate_limiter(struct cell *c, struct scheduler *s) {
- /* If this cell is already marked for drift, quit early. */
+ /* If this cell is already marked for limiting, quit early. */
if (c->hydro.do_limiter) return;
/* Mark this cell for limiting. */
@@ -1990,6 +2062,7 @@ void cell_activate_hydro_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_hydro_sorts_up(c, s);
}
@@ -2000,7 +2073,7 @@ void cell_activate_hydro_sorts(struct cell *c, int sid, struct scheduler *s) {
*/
void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
- if (c == c->super) {
+ if (c == c->hydro.super) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->stars.sorts == NULL)
error("Trying to activate un-existing c->stars.sorts");
@@ -2008,7 +2081,6 @@ void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
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 {
@@ -2017,16 +2089,13 @@ void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
parent != NULL && !parent->stars.do_sub_sort;
parent = parent->parent) {
parent->stars.do_sub_sort = 1;
- if (parent == c->super) {
+ if (parent == c->hydro.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);
- }
+ if (parent->nodeID == engine_rank) cell_activate_drift_spart(parent, s);
break;
}
}
@@ -2053,6 +2122,7 @@ void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s) {
/* 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);
}
@@ -2437,10 +2507,14 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
/* Store the current dx_max and h_max values. */
ci->stars.dx_max_part_old = ci->stars.dx_max_part;
ci->stars.h_max_old = ci->stars.h_max;
+ ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
+ ci->hydro.h_max_old = ci->hydro.h_max;
if (cj != NULL) {
cj->stars.dx_max_part_old = cj->stars.dx_max_part;
cj->stars.h_max_old = cj->stars.h_max;
+ cj->hydro.dx_max_part_old = cj->hydro.dx_max_part;
+ cj->hydro.h_max_old = cj->hydro.h_max;
}
/* Self interaction? */
@@ -2478,17 +2552,13 @@ 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;
- 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];
int sid = space_getsid(s->space, &ci, &cj, shift);
/* recurse? */
- if (cell_can_recurse_in_pair_stars_task(ci) &&
- cell_can_recurse_in_pair_stars_task(cj)) {
+ if (cell_can_recurse_in_pair_stars_task(ci, cj) &&
+ cell_can_recurse_in_pair_stars_task(cj, ci)) {
/* Different types of flags. */
switch (sid) {
@@ -2767,8 +2837,7 @@ 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) && cj->hydro.count != 0 &&
- ci->stars.count != 0) {
+ if (cell_is_active_stars(ci, e)) {
/* 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);
@@ -2785,8 +2854,7 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
cell_activate_stars_sorts(ci, sid, s);
}
- if (cell_is_active_stars(cj, e) && ci->hydro.count != 0 &&
- cj->stars.count != 0) {
+ if (cell_is_active_stars(cj, e)) {
/* 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);
@@ -3030,7 +3098,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
}
/* Store current values of dx_max and h_max. */
else if (t->type == task_type_sub_pair || t->type == task_type_sub_self) {
- cell_activate_subcell_hydro_tasks(t->ci, t->cj, s);
+ cell_activate_subcell_hydro_tasks(ci, cj, s);
}
}
@@ -3156,7 +3224,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
- if (c->end_force != NULL) scheduler_activate(s, c->end_force);
+ if (c->hydro.end_force != NULL) scheduler_activate(s, c->hydro.end_force);
if (c->hydro.cooling != NULL) scheduler_activate(s, c->hydro.cooling);
if (c->hydro.star_formation != NULL)
scheduler_activate(s, c->hydro.star_formation);
@@ -3304,11 +3372,11 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
- if (c->end_force != NULL) scheduler_activate(s, c->end_force);
if (c->grav.down != NULL) scheduler_activate(s, c->grav.down);
if (c->grav.down_in != NULL) scheduler_activate(s, c->grav.down_in);
if (c->grav.mesh != NULL) scheduler_activate(s, c->grav.mesh);
if (c->grav.long_range != NULL) scheduler_activate(s, c->grav.long_range);
+ if (c->grav.end_force != NULL) scheduler_activate(s, c->grav.end_force);
if (c->logger != NULL) scheduler_activate(s, c->logger);
/* Subgrid tasks */
@@ -3334,9 +3402,14 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
struct engine *e = s->space->e;
+ const int with_feedback = (e->policy & engine_policy_feedback);
const int nodeID = e->nodeID;
int rebuild = 0;
+ if (!with_feedback && c->stars.drift != NULL && cell_is_active_stars(c, e)) {
+ cell_activate_drift_spart(c, s);
+ }
+
/* Un-skip the density tasks involved with this cell. */
for (struct link *l = c->stars.density; l != NULL; l = l->next) {
struct task *t = l->t;
@@ -3344,60 +3417,69 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
struct cell *cj = t->cj;
const int ci_active = cell_is_active_stars(ci, e);
const int cj_active = (cj != NULL) ? cell_is_active_stars(cj, e) : 0;
+#ifdef WITH_MPI
+ const int ci_nodeID = ci->nodeID;
+ const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
+#else
+ const int ci_nodeID = nodeID;
+ const int cj_nodeID = nodeID;
+#endif
+
+ /* Activate the drifts */
+ if (t->type == task_type_self && ci_active) {
+ cell_activate_drift_part(ci, s);
+ cell_activate_drift_spart(ci, s);
+ }
/* Only activate tasks that involve a local active cell. */
- if ((ci_active && ci->nodeID == nodeID) ||
- (cj_active && cj->nodeID == nodeID)) {
+ if ((ci_active || cj_active) &&
+ (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
+
scheduler_activate(s, t);
- /* Activate drifts */
- if (t->type == task_type_self) {
- if (ci->nodeID == nodeID) {
- cell_activate_drift_part(ci, s);
- cell_activate_drift_spart(ci, s);
- }
- }
+ if (t->type == task_type_pair) {
- /* Set the correct sorting flags and activate hydro drifts */
- else if (t->type == task_type_pair) {
/* 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;
+ if (ci_active) {
+ /* 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);
- cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
+ /* hydro for cj */
+ atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
+ cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
- /* Activate the drift tasks. */
- if (ci->nodeID == nodeID) cell_activate_drift_spart(ci, s);
- if (cj->nodeID == nodeID) cell_activate_drift_part(cj, s);
+ /* Activate the drift tasks. */
+ 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_stars_sorts(ci, t->flags, s);
- cell_activate_hydro_sorts(cj, t->flags, s);
+ /* Check the sorts and activate them if needed. */
+ 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;
+ if (cj_active) {
+ /* hydro for ci */
+ atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
+ ci->hydro.dx_max_sort_old = ci->hydro.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);
+ /* stars for cj */
+ atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
+ cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
- /* Check the sorts and activate them if needed. */
- cell_activate_hydro_sorts(ci, t->flags, s);
- cell_activate_stars_sorts(cj, t->flags, s);
+ /* 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) {
- cell_activate_subcell_stars_tasks(t->ci, t->cj, s);
+ cell_activate_subcell_stars_tasks(ci, cj, s);
}
}
@@ -3407,99 +3489,102 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
/* Check whether there was too much particle motion, i.e. the
cell neighbour conditions were violated. */
if (cell_need_rebuild_for_stars_pair(ci, cj)) rebuild = 1;
+ if (cell_need_rebuild_for_stars_pair(cj, ci)) rebuild = 1;
#ifdef WITH_MPI
- error("MPI with stars not implemented");
- /* /\* Activate the send/recv tasks. *\/ */
- /* if (ci->nodeID != nodeID) { */
-
- /* /\* If the local cell is active, receive data from the foreign cell.
- * *\/ */
- /* if (cj_active) { */
- /* scheduler_activate(s, ci->hydro.recv_xv); */
- /* if (ci_active) { */
- /* scheduler_activate(s, ci->hydro.recv_rho); */
-
- /* } */
- /* } */
+ /* Activate the send/recv tasks. */
+ if (ci_nodeID != nodeID) {
- /* /\* If the foreign cell is active, we want its ti_end values. *\/ */
- /* if (ci_active) scheduler_activate(s, ci->mpi.recv_ti); */
+ if (cj_active) {
+ scheduler_activate(s, ci->mpi.hydro.recv_xv);
+ scheduler_activate(s, ci->mpi.hydro.recv_rho);
- /* /\* Is the foreign cell active and will need stuff from us? *\/ */
- /* if (ci_active) { */
+ /* If the local cell is active, more stuff will be needed. */
+ scheduler_activate_send(s, cj->mpi.stars.send, ci_nodeID);
+ cell_activate_drift_spart(cj, s);
- /* scheduler_activate_send(s, cj->hydro.send_xv, ci->nodeID); */
+ /* If the local cell is active, send its ti_end values. */
+ scheduler_activate_send(s, cj->mpi.send_ti, ci_nodeID);
+ }
- /* /\* Drift the cell which will be sent; note that not all sent */
- /* particles will be drifted, only those that are needed. *\/ */
- /* cell_activate_drift_part(cj, s); */
+ if (ci_active) {
+ scheduler_activate(s, ci->mpi.stars.recv);
- /* /\* If the local cell is also active, more stuff will be needed.
- * *\/ */
- /* if (cj_active) { */
- /* scheduler_activate_send(s, cj->hydro.send_rho, ci->nodeID); */
+ /* If the foreign cell is active, we want its ti_end values. */
+ scheduler_activate(s, ci->mpi.recv_ti);
- /* } */
- /* } */
+ /* Is the foreign cell active and will need stuff from us? */
+ scheduler_activate_send(s, cj->mpi.hydro.send_xv, ci_nodeID);
+ scheduler_activate_send(s, cj->mpi.hydro.send_rho, ci_nodeID);
- /* /\* If the local cell is active, send its ti_end values. *\/ */
- /* if (cj_active) scheduler_activate_send(s, cj->mpi.send_ti,
- * ci->nodeID);
- */
+ /* Drift the cell which will be sent; note that not all sent
+ particles will be drifted, only those that are needed. */
+ cell_activate_drift_part(cj, s);
+ }
- /* } else if (cj->nodeID != nodeID) { */
+ } else if (cj_nodeID != nodeID) {
- /* /\* If the local cell is active, receive data from the foreign cell.
- * *\/ */
- /* if (ci_active) { */
- /* scheduler_activate(s, cj->hydro.recv_xv); */
- /* if (cj_active) { */
- /* scheduler_activate(s, cj->hydro.recv_rho); */
+ /* If the local cell is active, receive data from the foreign cell. */
+ if (ci_active) {
+ scheduler_activate(s, cj->mpi.hydro.recv_xv);
+ scheduler_activate(s, cj->mpi.hydro.recv_rho);
- /* } */
- /* } */
+ /* If the local cell is active, more stuff will be needed. */
+ scheduler_activate_send(s, ci->mpi.stars.send, cj_nodeID);
+ cell_activate_drift_spart(ci, s);
- /* /\* If the foreign cell is active, we want its ti_end values. *\/ */
- /* if (cj_active) scheduler_activate(s, cj->mpi.recv_ti); */
+ /* If the local cell is active, send its ti_end values. */
+ scheduler_activate_send(s, ci->mpi.send_ti, cj_nodeID);
+ }
- /* /\* Is the foreign cell active and will need stuff from us? *\/ */
- /* if (cj_active) { */
+ if (cj_active) {
+ scheduler_activate(s, cj->mpi.stars.recv);
- /* scheduler_activate_send(s, ci->hydro.send_xv, cj->nodeID); */
+ /* If the foreign cell is active, we want its ti_end values. */
+ scheduler_activate(s, cj->mpi.recv_ti);
- /* /\* Drift the cell which will be sent; note that not all sent */
- /* particles will be drifted, only those that are needed. *\/ */
- /* cell_activate_drift_part(ci, s); */
+ /* Is the foreign cell active and will need stuff from us? */
+ scheduler_activate_send(s, ci->mpi.hydro.send_xv, cj_nodeID);
+ scheduler_activate_send(s, ci->mpi.hydro.send_rho, cj_nodeID);
- /* /\* If the local cell is also active, more stuff will be needed.
- * *\/ */
- /* if (ci_active) { */
+ /* Drift the cell which will be sent; note that not all sent
+ particles will be drifted, only those that are needed. */
+ cell_activate_drift_part(ci, s);
+ }
+ }
+#endif
+ }
+ }
- /* scheduler_activate_send(s, ci->hydro.send_rho, cj->nodeID); */
+ /* Un-skip the feedback tasks involved with this cell. */
+ for (struct link *l = c->stars.feedback; l != NULL; l = l->next) {
+ struct task *t = l->t;
+ struct cell *ci = t->ci;
+ struct cell *cj = t->cj;
+ const int ci_active = cell_is_active_stars(ci, e);
+ const int cj_active = (cj != NULL) ? cell_is_active_stars(cj, e) : 0;
+#ifdef WITH_MPI
+ const int ci_nodeID = ci->nodeID;
+ const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
+#else
+ const int ci_nodeID = nodeID;
+ const int cj_nodeID = nodeID;
+#endif
- /* } */
- /* } */
+ if ((ci_active && cj_nodeID == nodeID) ||
+ (cj_active && ci_nodeID == nodeID)) {
+ scheduler_activate(s, t);
- /* /\* If the local cell is active, send its ti_end values. *\/ */
- /* if (ci_active) scheduler_activate_send(s, ci->mpi.send_ti,
- * cj->nodeID);
- */
- /* } */
-#endif
+ /* Nothing more to do here, all drifts and sorts activated above */
}
}
/* Unskip all the other task types. */
if (c->nodeID == nodeID && cell_is_active_stars(c, e)) {
- /* Un-skip the feedback tasks involved with this cell. */
- for (struct link *l = c->stars.feedback; l != NULL; l = l->next)
- scheduler_activate(s, l->t);
-
- if (c->stars.ghost_in != NULL) scheduler_activate(s, c->stars.ghost_in);
- if (c->stars.ghost_out != NULL) scheduler_activate(s, c->stars.ghost_out);
if (c->stars.ghost != NULL) scheduler_activate(s, c->stars.ghost);
+ if (c->stars.stars_in != NULL) scheduler_activate(s, c->stars.stars_in);
+ if (c->stars.stars_out != NULL) scheduler_activate(s, c->stars.stars_out);
if (c->logger != NULL) scheduler_activate(s, c->logger);
}
@@ -3514,8 +3599,9 @@ int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s) {
*/
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 || c->grav.nr_mm_tasks > 0)) super = c;
+ /* Are we in a cell which is either the hydro or gravity super? */
+ if (super == NULL && (c->hydro.super != NULL || c->grav.super != NULL))
+ super = c;
/* Set the super-cell */
c->super = super;
@@ -3832,15 +3918,9 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
const int periodic = e->s->periodic;
const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
const int with_cosmology = (e->policy & engine_policy_cosmology);
- const float stars_h_max = e->stars_properties->h_max;
const integertime_t ti_old_gpart = c->grav.ti_old_part;
const integertime_t ti_current = e->ti_current;
struct gpart *const gparts = c->grav.parts;
- struct spart *const sparts = c->stars.parts;
-
- float dx_max = 0.f, dx2_max = 0.f;
- float dx_max_sort = 0.0f, dx2_max_sort = 0.f;
- float cell_h_max = 0.f;
/* Drift irrespective of cell flags? */
force |= c->grav.do_drift;
@@ -3878,19 +3958,9 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
/* Recurse */
cell_drift_gpart(cp, e, force);
-
- /* Update */
- dx_max = max(dx_max, cp->stars.dx_max_part);
- dx_max_sort = max(dx_max_sort, cp->stars.dx_max_sort);
- cell_h_max = max(cell_h_max, cp->stars.h_max);
}
}
- /* Store the values */
- c->stars.h_max = cell_h_max;
- c->stars.dx_max_part = dx_max;
- c->stars.dx_max_sort = dx_max_sort;
-
/* Update the time of the last drift */
c->grav.ti_old_part = ti_current;
@@ -3948,6 +4018,100 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
}
}
+ /* Update the time of the last drift */
+ c->grav.ti_old_part = ti_current;
+ }
+
+ /* Clear the drift flags. */
+ c->grav.do_drift = 0;
+ c->grav.do_sub_drift = 0;
+}
+
+/**
+ * @brief Recursively drifts the #spart in a cell hierarchy.
+ *
+ * @param c The #cell.
+ * @param e The #engine (to get ti_current).
+ * @param force Drift the particles irrespective of the #cell flags.
+ */
+void cell_drift_spart(struct cell *c, const struct engine *e, int force) {
+
+ const int periodic = e->s->periodic;
+ const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const float stars_h_max = e->hydro_properties->h_max;
+ const float stars_h_min = e->hydro_properties->h_min;
+ const integertime_t ti_old_spart = c->stars.ti_old_part;
+ const integertime_t ti_current = e->ti_current;
+ struct spart *const sparts = c->stars.parts;
+
+ float dx_max = 0.f, dx2_max = 0.f;
+ float dx_max_sort = 0.0f, dx2_max_sort = 0.f;
+ float cell_h_max = 0.f;
+
+ /* Drift irrespective of cell flags? */
+ force |= c->stars.do_drift;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that we only drift local cells. */
+ if (c->nodeID != engine_rank) error("Drifting a foreign cell is nope.");
+
+ /* Check that we are actually going to move forward. */
+ if (ti_current < ti_old_spart) error("Attempt to drift to the past");
+#endif
+
+ /* Early abort? */
+ if (c->stars.count == 0) {
+
+ /* Clear the drift flags. */
+ c->stars.do_drift = 0;
+ c->stars.do_sub_drift = 0;
+
+ /* Update the time of the last drift */
+ c->stars.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->stars.do_sub_drift)) {
+
+ /* Loop over the progeny and collect their data. */
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+ struct cell *cp = c->progeny[k];
+
+ /* Recurse */
+ cell_drift_spart(cp, e, force);
+
+ /* Update */
+ dx_max = max(dx_max, cp->stars.dx_max_part);
+ dx_max_sort = max(dx_max_sort, cp->stars.dx_max_sort);
+ cell_h_max = max(cell_h_max, cp->stars.h_max);
+ }
+ }
+
+ /* Store the values */
+ c->stars.h_max = cell_h_max;
+ c->stars.dx_max_part = dx_max;
+ c->stars.dx_max_sort = dx_max_sort;
+
+ /* Update the time of the last drift */
+ c->stars.ti_old_part = ti_current;
+
+ } else if (!c->split && force && ti_current > ti_old_spart) {
+
+ /* Drift from the last time the cell was drifted to the current time */
+ double dt_drift;
+ if (with_cosmology) {
+ dt_drift =
+ cosmology_get_drift_factor(e->cosmology, ti_old_spart, ti_current);
+ } else {
+ dt_drift = (ti_current - ti_old_spart) * e->time_base;
+ }
+
/* Loop over all the star particles in the cell */
const size_t nr_sparts = c->stars.count;
for (size_t k = 0; k < nr_sparts; k++) {
@@ -3959,7 +4123,7 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
if (spart_is_inhibited(sp, e)) continue;
/* Drift... */
- drift_spart(sp, dt_drift, ti_old_gpart, ti_current);
+ drift_spart(sp, dt_drift, ti_old_spart, ti_current);
#ifdef SWIFT_DEBUG_CHECKS
/* Make sure the particle does not drift by more than a box length. */
@@ -3991,6 +4155,7 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
/* Limit h to within the allowed range */
sp->h = min(sp->h, stars_h_max);
+ sp->h = max(sp->h, stars_h_min);
/* Compute (square of) motion since last cell construction */
const float dx2 = sp->x_diff[0] * sp->x_diff[0] +
@@ -4007,7 +4172,11 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
/* Maximal smoothing length */
cell_h_max = max(cell_h_max, sp->h);
- } /* Note: no need to compute dx_max as all spart have a gpart */
+ /* Get ready for a density calculation */
+ if (spart_is_active(sp, e)) {
+ stars_init_spart(sp);
+ }
+ }
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
@@ -4019,12 +4188,12 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
c->stars.dx_max_sort = dx_max_sort;
/* Update the time of the last drift */
- c->grav.ti_old_part = ti_current;
+ c->stars.ti_old_part = ti_current;
}
/* Clear the drift flags. */
- c->grav.do_drift = 0;
- c->grav.do_sub_drift = 0;
+ c->stars.do_drift = 0;
+ c->stars.do_sub_drift = 0;
}
/**
@@ -4099,6 +4268,42 @@ void cell_drift_multipole(struct cell *c, const struct engine *e) {
c->grav.ti_old_multipole = ti_current;
}
+/**
+ * @brief Resets all the sorting properties for the stars in a given cell
+ * hierarchy.
+ *
+ * @param c The #cell to clean.
+ * @param is_super Is this a super-cell?
+ */
+void cell_clear_stars_sort_flags(struct cell *c, const int is_super) {
+
+ /* Recurse if possible */
+ if (c->split) {
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ cell_clear_stars_sort_flags(c->progeny[k], /*is_super=*/0);
+ }
+
+ /* Free the sorted array at the level where it was allocated */
+ if (is_super) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c != c->hydro.super) error("Cell is not a super-cell!!!");
+#endif
+
+ for (int i = 0; i < 13; i++) {
+ free(c->stars.sort[i]);
+ }
+ }
+
+ /* Indicate that the cell is not sorted and cancel the pointer sorting arrays.
+ */
+ c->stars.sorted = 0;
+ for (int i = 0; i < 13; i++) {
+ c->stars.sort[i] = NULL;
+ }
+}
+
/**
* @brief Recursively checks that all particles in a cell have a time-step
*/
diff --git a/src/cell.h b/src/cell.h
index c5fbc9b8c02b0e008d337189fdbf582faf4fa600..b13889abda54f77ec2030bf2262dca30440bdc4e 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -150,6 +150,12 @@ struct pcell {
/*! Minimal integer end-of-timestep in this cell for stars tasks */
integertime_t ti_end_min;
+ /*! Maximal integer end-of-timestep in this cell for stars tasks */
+ integertime_t ti_end_max;
+
+ /*! Integer time of the last drift of the #spart in this cell */
+ integertime_t ti_old_part;
+
} stars;
/*! Maximal depth in that part of the tree */
@@ -198,12 +204,15 @@ struct pcell_step {
/*! Stars variables */
struct {
- /*! Maximal distance any #part has travelled since last rebuild */
- float dx_max_part;
-
/*! Minimal integer end-of-timestep in this cell (stars) */
integertime_t ti_end_min;
+ /*! Maximal integer end-of-timestep in this cell (stars) */
+ integertime_t ti_end_max;
+
+ /*! Maximal distance any #part has travelled since last rebuild */
+ float dx_max_part;
+
} stars;
};
@@ -278,6 +287,9 @@ struct cell {
/*! The extra ghost task for complex hydro schemes */
struct task *extra_ghost;
+ /*! The task to end the force calculation */
+ struct task *end_force;
+
/*! Task for cooling */
struct task *cooling;
@@ -409,6 +421,9 @@ struct cell {
/*! Task propagating the multipole to the particles */
struct task *down;
+ /*! The task to end the force calculation */
+ struct task *end_force;
+
/*! Minimum end of (integer) time step in this cell for gravity tasks. */
integertime_t ti_end_min;
@@ -466,12 +481,6 @@ struct cell {
/*! Pointer to the #spart data. */
struct spart *parts;
- /*! Dependency implicit task for the star ghost (in->ghost->out)*/
- struct task *ghost_in;
-
- /*! Dependency implicit task for the star ghost (in->ghost->out)*/
- struct task *ghost_out;
-
/*! The star ghost task itself */
struct task *ghost;
@@ -481,12 +490,25 @@ struct cell {
/*! Linked list of the tasks computing this cell's star feedback. */
struct link *feedback;
- /*! The task computing this cell's sorts. */
+ /*! The task computing this cell's sorts before the density. */
struct task *sorts;
+ /*! The drift task for sparts */
+ struct task *drift;
+
+ /*! Implicit tasks marking the entry of the stellar physics block of tasks
+ */
+ struct task *stars_in;
+
+ /*! Implicit tasks marking the exit of the stellar physics block of tasks */
+ struct task *stars_out;
+
/*! Max smoothing length in this cell. */
double h_max;
+ /*! Last (integer) time the cell's spart were drifted forward in time. */
+ integertime_t ti_old_part;
+
/*! Spin lock for various uses (#spart case). */
swift_lock_type lock;
@@ -529,6 +551,13 @@ struct cell {
/*! Maximum end of (integer) time step in this cell for gravity tasks. */
integertime_t ti_end_min;
+ /*! Maximum end of (integer) time step in this cell for star tasks. */
+ integertime_t ti_end_max;
+
+ /*! Maximum beginning of (integer) time step in this cell for star tasks.
+ */
+ integertime_t ti_beg_max;
+
/*! Number of #spart updated in this cell. */
int updated;
@@ -538,6 +567,12 @@ struct cell {
/*! Is the #spart data of this cell being used in a sub-cell? */
int hold;
+ /*! Does this cell need to be drifted (stars)? */
+ char do_drift;
+
+ /*! Do any of this cell's sub-cells need to be drifted (stars)? */
+ char do_sub_drift;
+
#ifdef SWIFT_DEBUG_CHECKS
/*! Last (integer) time the cell's sort arrays were updated. */
integertime_t ti_sort;
@@ -580,11 +615,18 @@ struct cell {
} grav;
struct {
+ /* Task receiving spart data. */
+ struct task *recv;
- /* Task receiving gpart data. */
+ /* Linked list for sending spart data. */
+ struct link *send;
+ } stars;
+
+ struct {
+ /* Task receiving limiter data. */
struct task *recv;
- /* Linked list for sending gpart data. */
+ /* Linked list for sending limiter data. */
struct link *send;
} limiter;
@@ -609,9 +651,6 @@ struct cell {
} mpi;
#endif
- /*! The task to end the force calculation */
- struct task *end_force;
-
/*! The first kick task */
struct task *kick1;
@@ -702,6 +741,7 @@ void cell_check_foreign_multipole(const struct cell *c);
void cell_clean(struct cell *c);
void cell_check_part_drift_point(struct cell *c, void *data);
void cell_check_gpart_drift_point(struct cell *c, void *data);
+void cell_check_spart_drift_point(struct cell *c, void *data);
void cell_check_multipole_drift_point(struct cell *c, void *data);
void cell_reset_task_counters(struct cell *c);
int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s);
@@ -710,6 +750,7 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
void cell_drift_part(struct cell *c, const struct engine *e, int force);
void cell_drift_gpart(struct cell *c, const struct engine *e, int force);
+void cell_drift_spart(struct cell *c, const struct engine *e, int force);
void cell_drift_multipole(struct cell *c, const struct engine *e);
void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
@@ -733,6 +774,7 @@ void cell_clear_limiter_flags(struct cell *c, void *data);
void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data);
void cell_check_spart_pos(const struct cell *c,
const struct spart *global_sparts);
+void cell_clear_stars_sort_flags(struct cell *c, const int is_super);
int cell_has_tasks(struct cell *c);
void cell_remove_part(const struct engine *e, struct cell *c, struct part *p,
struct xpart *xp);
@@ -854,17 +896,22 @@ cell_can_recurse_in_self_hydro_task(const struct cell *c) {
* @brief Can a sub-pair star task recurse to a lower level based
* on the status of the particles in the cell.
*
- * @param c The #cell.
+ * @param ci The #cell with stars.
+ * @param cj The #cell with hydro parts.
*/
__attribute__((always_inline)) INLINE static int
-cell_can_recurse_in_pair_stars_task(const struct cell *c) {
+cell_can_recurse_in_pair_stars_task(const struct cell *ci,
+ const struct cell *cj) {
/* Is the cell split ? */
/* If so, is the cut-off radius plus the max distance the parts have moved */
/* smaller than the sub-cell sizes ? */
/* Note: We use the _old values as these might have been updated by a drift */
- return c->split && ((kernel_gamma * c->stars.h_max_old +
- c->stars.dx_max_part_old) < 0.5f * c->dmin);
+ return ci->split && cj->split &&
+ ((kernel_gamma * ci->stars.h_max_old + ci->stars.dx_max_part_old) <
+ 0.5f * ci->dmin) &&
+ ((kernel_gamma * cj->hydro.h_max_old + cj->hydro.dx_max_part_old) <
+ 0.5f * cj->dmin);
}
/**
@@ -877,7 +924,8 @@ __attribute__((always_inline)) INLINE static int
cell_can_recurse_in_self_stars_task(const struct cell *c) {
/* Is the cell split and not smaller than the smoothing length? */
- return c->split && (kernel_gamma * c->stars.h_max_old < 0.5f * c->dmin);
+ return c->split && (kernel_gamma * c->stars.h_max_old < 0.5f * c->dmin) &&
+ (kernel_gamma * c->hydro.h_max_old < 0.5f * c->dmin);
}
/**
@@ -895,7 +943,8 @@ __attribute__((always_inline)) INLINE static int cell_can_split_pair_hydro_task(
/* Note that since tasks are create after a rebuild no need to take */
/* into account any part motion (i.e. dx_max == 0 here) */
return c->split &&
- (space_stretch * kernel_gamma * c->hydro.h_max < 0.5f * c->dmin);
+ (space_stretch * kernel_gamma * c->hydro.h_max < 0.5f * c->dmin) &&
+ (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin);
}
/**
@@ -913,42 +962,7 @@ __attribute__((always_inline)) INLINE static int cell_can_split_self_hydro_task(
/* Note: No need for more checks here as all the sub-pairs and sub-self */
/* tasks will be created. So no need to check for h_max */
return c->split &&
- (space_stretch * kernel_gamma * c->hydro.h_max < 0.5f * c->dmin);
-}
-
-/**
- * @brief Can a pair stars task associated with a cell be split into smaller
- * sub-tasks.
- *
- * @param c The #cell.
- */
-__attribute__((always_inline)) INLINE static int cell_can_split_pair_stars_task(
- const struct cell *c) {
-
- /* Is the cell split ? */
- /* If so, is the cut-off radius with some leeway smaller than */
- /* the sub-cell sizes ? */
- /* Note that since tasks are create after a rebuild no need to take */
- /* into account any part motion (i.e. dx_max == 0 here) */
- return c->split &&
- (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin);
-}
-
-/**
- * @brief Can a self stars task associated with a cell be split into smaller
- * sub-tasks.
- *
- * @param c The #cell.
- */
-__attribute__((always_inline)) INLINE static int cell_can_split_self_stars_task(
- const struct cell *c) {
-
- /* Is the cell split ? */
- /* If so, is the cut-off radius with some leeway smaller than */
- /* the sub-cell sizes ? */
- /* Note: No need for more checks here as all the sub-pairs and sub-self */
- /* tasks will be created. So no need to check for h_max */
- return c->split &&
+ (space_stretch * kernel_gamma * c->hydro.h_max < 0.5f * c->dmin) &&
(space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin);
}
@@ -992,15 +1006,16 @@ cell_need_rebuild_for_hydro_pair(const struct cell *ci, const struct cell *cj) {
/* Is the cut-off radius plus the max distance the parts in both cells have */
/* moved larger than the cell size ? */
/* Note ci->dmin == cj->dmin */
- return (kernel_gamma * max(ci->hydro.h_max, cj->hydro.h_max) +
- ci->hydro.dx_max_part + cj->hydro.dx_max_part >
- cj->dmin);
+ if (kernel_gamma * max(ci->hydro.h_max, cj->hydro.h_max) +
+ ci->hydro.dx_max_part + cj->hydro.dx_max_part >
+ cj->dmin) {
+ return 1;
+ }
+ return 0;
}
-
/**
* @brief Have star particles in a pair of cells moved too much and require a
- * rebuild
- * ?
+ * rebuild?
*
* @param ci The first #cell.
* @param cj The second #cell.
@@ -1011,9 +1026,12 @@ cell_need_rebuild_for_stars_pair(const struct cell *ci, const struct cell *cj) {
/* Is the cut-off radius plus the max distance the parts in both cells have */
/* moved larger than the cell size ? */
/* Note ci->dmin == cj->dmin */
- return (kernel_gamma * max(ci->stars.h_max, cj->stars.h_max) +
- ci->stars.dx_max_part + cj->stars.dx_max_part >
- cj->dmin);
+ if (kernel_gamma * max(ci->stars.h_max, cj->hydro.h_max) +
+ ci->stars.dx_max_part + cj->hydro.dx_max_part >
+ cj->dmin) {
+ return 1;
+ }
+ return 0;
}
/**
diff --git a/src/debug.c b/src/debug.c
index 2c1d81b676520c4823a27f3a7a8cc617585ca5f2..d2aff378a174ade46b62a3931f78394a0f41ca41 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -181,6 +181,14 @@ int checkSpacehmax(struct space *s) {
}
}
+ float cell_stars_h_max = 0.0f;
+ for (int k = 0; k < s->nr_cells; k++) {
+ if (s->cells_top[k].nodeID == s->e->nodeID &&
+ s->cells_top[k].stars.h_max > cell_stars_h_max) {
+ cell_stars_h_max = s->cells_top[k].stars.h_max;
+ }
+ }
+
/* Now all particles. */
float part_h_max = 0.0f;
for (size_t k = 0; k < s->nr_parts; k++) {
@@ -189,10 +197,21 @@ int checkSpacehmax(struct space *s) {
}
}
+ /* Now all the sparticles. */
+ float spart_h_max = 0.0f;
+ for (size_t k = 0; k < s->nr_sparts; k++) {
+ if (s->sparts[k].h > spart_h_max) {
+ spart_h_max = s->sparts[k].h;
+ }
+ }
+
/* If within some epsilon we are OK. */
- if (fabsf(cell_h_max - part_h_max) <= FLT_EPSILON) return 1;
+ if (fabsf(cell_h_max - part_h_max) <= FLT_EPSILON &&
+ fabsf(cell_stars_h_max - spart_h_max) <= FLT_EPSILON)
+ return 1;
/* There is a problem. Hunt it down. */
+ /* part */
for (int k = 0; k < s->nr_cells; k++) {
if (s->cells_top[k].nodeID == s->e->nodeID) {
if (s->cells_top[k].hydro.h_max > part_h_max) {
@@ -209,6 +228,23 @@ int checkSpacehmax(struct space *s) {
}
}
+ /* spart */
+ for (int k = 0; k < s->nr_cells; k++) {
+ if (s->cells_top[k].nodeID == s->e->nodeID) {
+ if (s->cells_top[k].stars.h_max > spart_h_max) {
+ message("cell %d is inconsistent (%f > %f)", k,
+ s->cells_top[k].stars.h_max, spart_h_max);
+ }
+ }
+ }
+
+ for (size_t k = 0; k < s->nr_sparts; k++) {
+ if (s->sparts[k].h > cell_stars_h_max) {
+ message("spart %lld is inconsistent (%f > %f)", s->sparts[k].id,
+ s->sparts[k].h, cell_stars_h_max);
+ }
+ }
+
return 0;
}
@@ -227,6 +263,8 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
float h_max = 0.0f;
float dx_max = 0.0f;
+ float stars_h_max = 0.0f;
+ float stars_dx_max = 0.0f;
int result = 1;
const double loc_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
@@ -262,6 +300,33 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
dx_max = max(dx_max, sqrt(dx2));
}
+ const size_t nr_sparts = c->stars.count;
+ struct spart *sparts = c->stars.parts;
+ for (size_t k = 0; k < nr_sparts; k++) {
+
+ struct spart *const sp = &sparts[k];
+
+ if (sp->x[0] < loc_min[0] || sp->x[0] >= loc_max[0] ||
+ sp->x[1] < loc_min[1] || sp->x[1] >= loc_max[1] ||
+ sp->x[2] < loc_min[2] || sp->x[2] >= loc_max[2]) {
+
+ message(
+ "Inconsistent part position p->x=[%e %e %e], c->loc=[%e %e %e] "
+ "c->width=[%e %e %e]",
+ sp->x[0], sp->x[1], sp->x[2], c->loc[0], c->loc[1], c->loc[2],
+ c->width[0], c->width[1], c->width[2]);
+
+ result = 0;
+ }
+
+ const float dx2 = sp->x_diff[0] * sp->x_diff[0] +
+ sp->x_diff[1] * sp->x_diff[1] +
+ sp->x_diff[2] * sp->x_diff[2];
+
+ stars_h_max = max(stars_h_max, sp->h);
+ stars_dx_max = max(stars_dx_max, sqrt(dx2));
+ }
+
if (c->split) {
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
@@ -285,6 +350,19 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
result = 0;
}
+ if (c->stars.h_max != stars_h_max) {
+ message("%d Inconsistent stars_h_max: cell %f != parts %f", *depth,
+ c->stars.h_max, stars_h_max);
+ message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+ result = 0;
+ }
+ if (c->stars.dx_max_part != stars_dx_max) {
+ message("%d Inconsistent stars_dx_max: %f != %f", *depth,
+ c->stars.dx_max_part, stars_dx_max);
+ message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+ result = 0;
+ }
+
return result;
}
diff --git a/src/drift.h b/src/drift.h
index a4bdf9be74aade4fe0f1349544cf472363c81c99..7e874fe0ceabe5b091cc7c5bb53adbef2c9a3efd 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -28,6 +28,7 @@
#include "dimension.h"
#include "hydro.h"
#include "part.h"
+#include "stars.h"
/**
* @brief Perform the 'drift' operation on a #gpart.
@@ -138,6 +139,9 @@ __attribute__((always_inline)) INLINE static void drift_spart(
sp->x[1] += sp->v[1] * dt_drift;
sp->x[2] += sp->v[2] * dt_drift;
+ /* Predict the values of the extra fields */
+ stars_predict_extra(sp, dt_drift);
+
/* Compute offsets since last cell construction */
for (int k = 0; k < 3; k++) {
const float dx = sp->v[k] * dt_drift;
diff --git a/src/engine.c b/src/engine.c
index bc243bac7bee9d862071c8487134e2a3b7902902..c8fe9e9c777e494540a2b6a9f424157eba06e597 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2710,16 +2710,21 @@ void engine_skip_force_and_kick(struct engine *e) {
/* Skip everything that updates the particles */
if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
- t->type == task_type_kick1 || t->type == task_type_kick2 ||
- t->type == task_type_timestep ||
+ t->type == task_type_drift_spart || t->type == task_type_kick1 ||
+ t->type == task_type_kick2 || t->type == task_type_timestep ||
t->type == task_type_timestep_limiter ||
t->subtype == task_subtype_force ||
t->subtype == task_subtype_limiter || t->subtype == task_subtype_grav ||
- t->type == task_type_end_force ||
+ t->type == task_type_end_hydro_force ||
+ t->type == task_type_end_grav_force ||
t->type == task_type_grav_long_range || t->type == task_type_grav_mm ||
- t->type == task_type_grav_down || t->type == task_type_cooling ||
+ t->type == task_type_grav_down || t->type == task_type_grav_down_in ||
+ t->type == task_type_drift_gpart_out || t->type == task_type_cooling ||
+ t->type == task_type_stars_in || t->type == task_type_stars_out ||
t->type == task_type_star_formation ||
- t->type == task_type_extra_ghost || t->subtype == task_subtype_gradient)
+ t->type == task_type_extra_ghost ||
+ t->subtype == task_subtype_gradient ||
+ t->subtype == task_subtype_stars_feedback)
t->skip = 1;
}
@@ -2743,7 +2748,8 @@ void engine_skip_drift(struct engine *e) {
struct task *t = &tasks[i];
/* Skip everything that moves the particles */
- if (t->type == task_type_drift_part || t->type == task_type_drift_gpart)
+ if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
+ t->type == task_type_drift_spart)
t->skip = 1;
}
@@ -2757,7 +2763,6 @@ void engine_skip_drift(struct engine *e) {
* @param e The #engine.
*/
void engine_launch(struct engine *e) {
-
const ticks tic = getticks();
#ifdef SWIFT_DEBUG_CHECKS
@@ -3949,6 +3954,78 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
#endif
}
+#ifdef DEBUG_INTERACTIONS_STARS
+/**
+ * @brief Exchange the feedback counters between stars
+ * @param e The #engine.
+ */
+void engine_collect_stars_counter(struct engine *e) {
+
+#ifdef WITH_MPI
+ if (e->total_nr_sparts > 1e5) {
+ message("WARNING: too many sparts, skipping exchange of counters");
+ return;
+ }
+
+ /* Get number of sparticles for each rank */
+ size_t *n_sparts = (size_t *)malloc(e->nr_nodes * sizeof(size_t));
+
+ int err = MPI_Allgather(&e->s->nr_sparts_foreign, 1, MPI_UNSIGNED_LONG,
+ n_sparts, 1, MPI_UNSIGNED_LONG, MPI_COMM_WORLD);
+ if (err != MPI_SUCCESS) error("Communication failed");
+
+ /* Compute derivated quantities */
+ int total = 0;
+ int *n_sparts_int = (int *)malloc(e->nr_nodes * sizeof(int));
+ int *displs = (int *)malloc(e->nr_nodes * sizeof(int));
+ for (int i = 0; i < e->nr_nodes; i++) {
+ displs[i] = total;
+ total += n_sparts[i];
+ n_sparts_int[i] = n_sparts[i];
+ }
+
+ /* Get all sparticles */
+ struct spart *sparts = (struct spart *)malloc(total * sizeof(struct spart));
+ err = MPI_Allgatherv(e->s->sparts_foreign, e->s->nr_sparts_foreign,
+ spart_mpi_type, sparts, n_sparts_int, displs,
+ spart_mpi_type, MPI_COMM_WORLD);
+ if (err != MPI_SUCCESS) error("Communication failed");
+
+ /* Reset counters */
+ for (size_t i = 0; i < e->s->nr_sparts_foreign; i++) {
+ e->s->sparts_foreign[i].num_ngb_force = 0;
+ }
+
+ /* Update counters */
+ struct spart *local_sparts = e->s->sparts;
+ for (size_t i = 0; i < e->s->nr_sparts; i++) {
+ const long long id_i = local_sparts[i].id;
+
+ for (int j = 0; j < total; j++) {
+ const long long id_j = sparts[j].id;
+
+ if (id_j == id_i) {
+ if (j >= displs[engine_rank] &&
+ j < displs[engine_rank] + n_sparts_int[engine_rank]) {
+ error(
+ "Found a local spart in foreign cell ID=%lli: j=%i, displs=%i, "
+ "n_sparts=%i",
+ id_j, j, displs[engine_rank], n_sparts_int[engine_rank]);
+ }
+
+ local_sparts[i].num_ngb_force += sparts[j].num_ngb_force;
+ }
+ }
+ }
+
+ free(n_sparts);
+ free(n_sparts_in);
+ free(sparts);
+#endif
+}
+
+#endif
+
/**
* @brief Writes a snapshot with the current state of the engine
*
@@ -3985,6 +4062,10 @@ void engine_dump_snapshot(struct engine *e) {
}
#endif
+#ifdef DEBUG_INTERACTIONS_STARS
+ engine_collect_stars_counter(e);
+#endif
+
/* Dump... */
#if defined(HAVE_HDF5)
#if defined(WITH_MPI)
diff --git a/src/engine_drift.c b/src/engine_drift.c
index 7a842068b57813575c33dd670172059abb1e8fc0..1b0711619d68da02753f307190ca3a0624feecce 100644
--- a/src/engine_drift.c
+++ b/src/engine_drift.c
@@ -124,6 +124,54 @@ void engine_do_drift_all_gpart_mapper(void *map_data, int num_elements,
}
}
+/**
+ * @brief Mapper function to drift *all* the #spart to the current time.
+ *
+ * @param map_data An array of #cell%s.
+ * @param num_elements Chunk size.
+ * @param extra_data Pointer to an #engine.
+ */
+void engine_do_drift_all_spart_mapper(void *map_data, int num_elements,
+ void *extra_data) {
+
+ const struct engine *e = (const struct engine *)extra_data;
+ const int restarting = e->restarting;
+ struct space *s = e->s;
+ struct cell *cells_top;
+ int *local_cells_top;
+
+ if (restarting) {
+
+ /* When restarting, we loop over all top-level cells */
+ cells_top = (struct cell *)map_data;
+ local_cells_top = NULL;
+
+ } else {
+
+ /* In any other case, we use the list of local cells with tasks */
+ cells_top = s->cells_top;
+ local_cells_top = (int *)map_data;
+ }
+
+ for (int ind = 0; ind < num_elements; ind++) {
+
+ struct cell *c;
+
+ /* When restarting, the list of local cells does not
+ yet exist. We use the raw list of top-level cells instead */
+ if (restarting)
+ c = &cells_top[ind];
+ else
+ c = &cells_top[local_cells_top[ind]];
+
+ if (c->nodeID == e->nodeID) {
+
+ /* Drift all the particles */
+ cell_drift_spart(c, e, /* force the drift=*/1);
+ }
+ }
+}
+
/**
* @brief Mapper function to drift *all* the multipoles to the current time.
*
@@ -204,6 +252,11 @@ void engine_drift_all(struct engine *e, const int drift_mpoles) {
e->s->local_cells_top, e->s->nr_local_cells, sizeof(int),
/* default chunk */ 0, e);
}
+ if (e->s->nr_sparts > 0) {
+ threadpool_map(&e->threadpool, engine_do_drift_all_spart_mapper,
+ e->s->local_cells_top, e->s->nr_local_cells, sizeof(int),
+ /* default chunk */ 0, e);
+ }
if (drift_mpoles && (e->policy & engine_policy_self_gravity)) {
threadpool_map(&e->threadpool, engine_do_drift_all_multipole_mapper,
e->s->local_cells_with_tasks_top,
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index d1d3be6c12ee47e22da42b41456c822394ed2ba4..6340cb17eb574822957aa4a7e6180cda9a14b5d7 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -146,6 +146,7 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
0, ci, cj);
t_rho = scheduler_addtask(s, task_type_send, task_subtype_rho,
ci->mpi.tag, 0, ci, cj);
+
#ifdef EXTRA_HYDRO_LOOP
t_gradient = scheduler_addtask(s, task_type_send, task_subtype_gradient,
ci->mpi.tag, 0, ci, cj);
@@ -153,7 +154,7 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
#ifdef EXTRA_HYDRO_LOOP
- scheduler_addunlock(s, t_gradient, ci->super->kick2);
+ scheduler_addunlock(s, t_gradient, ci->hydro.super->hydro.end_force);
scheduler_addunlock(s, ci->hydro.super->hydro.extra_ghost, t_gradient);
@@ -169,7 +170,7 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
#else
/* The send_rho task should unlock the super_hydro-cell's kick task. */
- scheduler_addunlock(s, t_rho, ci->super->end_force);
+ scheduler_addunlock(s, t_rho, ci->hydro.super->hydro.end_force);
/* The send_rho task depends on the cell's ghost task. */
scheduler_addunlock(s, ci->hydro.super->hydro.ghost_out, t_rho);
@@ -179,6 +180,8 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
#endif
+ scheduler_addunlock(s, ci->hydro.super->hydro.drift, t_rho);
+
/* Drift before you send */
scheduler_addunlock(s, ci->hydro.super->hydro.drift, t_xv);
}
@@ -203,6 +206,65 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
#endif
}
+/**
+ * @brief Add send tasks for the stars pairs to a hierarchy of cells.
+ *
+ * @param e The #engine.
+ * @param ci The sending #cell.
+ * @param cj Dummy cell containing the nodeID of the receiving node.
+ * @param t_feedback The send_feed #task, if it has already been created.
+ */
+void engine_addtasks_send_stars(struct engine *e, struct cell *ci,
+ struct cell *cj, struct task *t_feedback) {
+
+#ifdef WITH_MPI
+
+ struct link *l = NULL;
+ struct scheduler *s = &e->sched;
+ const int nodeID = cj->nodeID;
+
+ /* Check if any of the density tasks are for the target node. */
+ for (l = ci->stars.density; l != NULL; l = l->next)
+ if (l->t->ci->nodeID == nodeID ||
+ (l->t->cj != NULL && l->t->cj->nodeID == nodeID))
+ break;
+
+ /* If so, attach send tasks. */
+ if (l != NULL) {
+
+ if (t_feedback == NULL) {
+
+ /* Make sure this cell is tagged. */
+ cell_ensure_tagged(ci);
+
+ /* Create the tasks and their dependencies? */
+ t_feedback = scheduler_addtask(s, task_type_send, task_subtype_spart,
+ ci->mpi.tag, 0, ci, cj);
+
+ /* The send_stars task should unlock the super_cell's kick task. */
+ scheduler_addunlock(s, t_feedback, ci->hydro.super->stars.stars_out);
+
+ /* Ghost before you send */
+ scheduler_addunlock(s, ci->hydro.super->stars.ghost, t_feedback);
+
+ /* Drift before you send */
+ scheduler_addunlock(s, ci->hydro.super->stars.drift, t_feedback);
+ }
+
+ engine_addlink(e, &ci->mpi.stars.send, t_feedback);
+ }
+
+ /* Recurse? */
+ if (ci->split)
+ for (int k = 0; k < 8; k++)
+ if (ci->progeny[k] != NULL)
+ engine_addtasks_send_stars(e, ci->progeny[k], cj, t_feedback);
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+#endif
+}
+
/**
* @brief Add send tasks for the time-step to a hierarchy of cells.
*
@@ -236,6 +298,12 @@ void engine_addtasks_send_timestep(struct engine *e, struct cell *ci,
(l->t->cj != NULL && l->t->cj->nodeID == nodeID))
break;
+ if (l == NULL)
+ for (l = ci->stars.density; l != NULL; l = l->next)
+ if (l->t->ci->nodeID == nodeID ||
+ (l->t->cj != NULL && l->t->cj->nodeID == nodeID))
+ break;
+
/* If found anything, attach send tasks. */
if (l != NULL) {
@@ -319,7 +387,10 @@ void engine_addtasks_recv_hydro(struct engine *e, struct cell *c,
c->mpi.hydro.recv_gradient = t_gradient;
/* Add dependencies. */
- if (c->hydro.sorts != NULL) scheduler_addunlock(s, t_xv, c->hydro.sorts);
+ if (c->hydro.sorts != NULL) {
+ scheduler_addunlock(s, t_xv, c->hydro.sorts);
+ scheduler_addunlock(s, c->hydro.sorts, t_rho);
+ }
for (struct link *l = c->hydro.density; l != NULL; l = l->next) {
scheduler_addunlock(s, t_xv, l->t);
@@ -330,13 +401,20 @@ void engine_addtasks_recv_hydro(struct engine *e, struct cell *c,
scheduler_addunlock(s, t_rho, l->t);
scheduler_addunlock(s, l->t, t_gradient);
}
- for (struct link *l = c->hydro.force; l != NULL; l = l->next)
+ for (struct link *l = c->hydro.force; l != NULL; l = l->next) {
scheduler_addunlock(s, t_gradient, l->t);
+ }
#else
- for (struct link *l = c->hydro.force; l != NULL; l = l->next)
+ for (struct link *l = c->hydro.force; l != NULL; l = l->next) {
scheduler_addunlock(s, t_rho, l->t);
+ }
#endif
+ /* Make sure the density has been computed before the stars compute theirs. */
+ for (struct link *l = c->stars.density; l != NULL; l = l->next) {
+ scheduler_addunlock(s, t_rho, l->t);
+ }
+
/* Recurse? */
if (c->split)
for (int k = 0; k < 8; k++)
@@ -348,6 +426,59 @@ void engine_addtasks_recv_hydro(struct engine *e, struct cell *c,
#endif
}
+/**
+ * @brief Add recv tasks for stars pairs to a hierarchy of cells.
+ *
+ * @param e The #engine.
+ * @param c The foreign #cell.
+ * @param t_feedback The recv_feed #task, if it has already been created.
+ */
+void engine_addtasks_recv_stars(struct engine *e, struct cell *c,
+ struct task *t_feedback) {
+
+#ifdef WITH_MPI
+ struct scheduler *s = &e->sched;
+
+ /* Have we reached a level where there are any stars tasks ? */
+ if (t_feedback == NULL && c->stars.density != NULL) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Make sure this cell has a valid tag. */
+ if (c->mpi.tag < 0) error("Trying to receive from untagged cell.");
+#endif // SWIFT_DEBUG_CHECKS
+
+ /* Create the tasks. */
+ t_feedback = scheduler_addtask(s, task_type_recv, task_subtype_spart,
+ c->mpi.tag, 0, c, NULL);
+ }
+
+ c->mpi.stars.recv = t_feedback;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->nodeID == e->nodeID) error("Local cell!");
+#endif
+ if (c->stars.sorts != NULL)
+ scheduler_addunlock(s, t_feedback, c->stars.sorts);
+
+ for (struct link *l = c->stars.density; l != NULL; l = l->next) {
+ scheduler_addunlock(s, l->t, t_feedback);
+ }
+
+ for (struct link *l = c->stars.feedback; l != NULL; l = l->next) {
+ scheduler_addunlock(s, t_feedback, l->t);
+ }
+
+ /* Recurse? */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ engine_addtasks_recv_stars(e, c->progeny[k], t_feedback);
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+#endif
+}
+
/**
* @brief Add recv tasks for gravity pairs to a hierarchy of cells.
*
@@ -446,6 +577,9 @@ void engine_addtasks_recv_timestep(struct engine *e, struct cell *c,
}
}
+ for (struct link *l = c->stars.feedback; l != NULL; l = l->next)
+ scheduler_addunlock(s, l->t, t_ti);
+
/* Recurse? */
if (c->split)
for (int k = 0; k < 8; k++)
@@ -473,8 +607,6 @@ void engine_addtasks_recv_timestep(struct engine *e, struct cell *c,
void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
struct scheduler *s = &e->sched;
- const int is_with_cooling = (e->policy & engine_policy_cooling);
- const int is_with_star_formation = (e->policy & engine_policy_star_formation);
const int with_limiter = (e->policy & engine_policy_limiter);
/* Are we in a super-cell ? */
@@ -499,35 +631,7 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
c->timestep = scheduler_addtask(s, task_type_timestep, task_subtype_none,
0, 0, c, NULL);
- /* Add the task finishing the force calculation */
- c->end_force = scheduler_addtask(s, task_type_end_force,
- task_subtype_none, 0, 0, c, NULL);
-
- /* Subgrid tasks */
- if (is_with_cooling && c->hydro.count_total > 0) {
-
- c->hydro.cooling = scheduler_addtask(s, task_type_cooling,
- task_subtype_none, 0, 0, c, NULL);
-
- scheduler_addunlock(s, c->end_force, c->hydro.cooling);
- scheduler_addunlock(s, c->hydro.cooling, c->kick2);
-
- } else {
- scheduler_addunlock(s, c->end_force, c->kick2);
- }
-
- if (is_with_star_formation && c->hydro.count_total > 0) {
-
- c->hydro.star_formation = scheduler_addtask(
- s, task_type_star_formation, task_subtype_none, 0, 0, c, NULL);
-
- scheduler_addunlock(s, c->kick2, c->hydro.star_formation);
- scheduler_addunlock(s, c->hydro.star_formation, c->timestep);
-
- } else {
- scheduler_addunlock(s, c->kick2, c->timestep);
- }
-
+ scheduler_addunlock(s, c->kick2, c->timestep);
scheduler_addunlock(s, c->timestep, c->kick1);
/* Time-step limiting */
@@ -581,6 +685,11 @@ void engine_make_hierarchical_tasks_gravity(struct engine *e, struct cell *c) {
c->grav.drift = scheduler_addtask(s, task_type_drift_gpart,
task_subtype_none, 0, 0, c, NULL);
+ c->grav.end_force = scheduler_addtask(s, task_type_end_grav_force,
+ task_subtype_none, 0, 0, c, NULL);
+
+ scheduler_addunlock(s, c->grav.end_force, c->super->kick2);
+
if (is_self_gravity) {
/* Initialisation of the multipoles */
@@ -612,7 +721,7 @@ void engine_make_hierarchical_tasks_gravity(struct engine *e, struct cell *c) {
if (periodic) scheduler_addunlock(s, c->grav.mesh, c->grav.down);
scheduler_addunlock(s, c->grav.init, c->grav.long_range);
scheduler_addunlock(s, c->grav.long_range, c->grav.down);
- scheduler_addunlock(s, c->grav.down, c->super->end_force);
+ scheduler_addunlock(s, c->grav.down, c->grav.super->grav.end_force);
/* Link in the implicit tasks */
scheduler_addunlock(s, c->grav.init, c->grav.init_out);
@@ -654,34 +763,6 @@ void engine_make_hierarchical_tasks_gravity(struct engine *e, struct cell *c) {
engine_make_hierarchical_tasks_gravity(e, c->progeny[k]);
}
-/**
- * @brief Recursively add non-implicit star ghost tasks to a cell hierarchy.
- */
-void engine_add_stars_ghosts(struct engine *e, struct cell *c,
- struct task *stars_ghost_in,
- struct task *stars_ghost_out) {
-
- /* Abort as there are no star particles here? */
- if (c->stars.count_total == 0) return;
-
- /* If we have reached the leaf OR have to few particles to play with*/
- if (!c->split || c->stars.count_total < engine_max_sparts_per_ghost) {
-
- /* Add the ghost task and its dependencies */
- struct scheduler *s = &e->sched;
- c->stars.ghost = scheduler_addtask(s, task_type_stars_ghost,
- task_subtype_none, 0, 0, c, NULL);
- scheduler_addunlock(s, stars_ghost_in, c->stars.ghost);
- scheduler_addunlock(s, c->stars.ghost, stars_ghost_out);
- } else {
- /* Keep recursing */
- for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL)
- engine_add_stars_ghosts(e, c->progeny[k], stars_ghost_in,
- stars_ghost_out);
- }
-}
-
/**
* @brief Recursively add non-implicit ghost tasks to a cell hierarchy.
*/
@@ -724,6 +805,10 @@ void engine_add_ghosts(struct engine *e, struct cell *c, struct task *ghost_in,
void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
struct scheduler *s = &e->sched;
+ const int with_stars = (e->policy & engine_policy_stars);
+ const int with_feedback = (e->policy & engine_policy_feedback);
+ const int with_cooling = (e->policy & engine_policy_cooling);
+ const int with_star_formation = (e->policy & engine_policy_star_formation);
/* Are we in a super-cell ? */
if (c->hydro.super == c) {
@@ -732,6 +817,11 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
c->hydro.sorts =
scheduler_addtask(s, task_type_sort, task_subtype_none, 0, 0, c, NULL);
+ if (with_feedback) {
+ c->stars.sorts = scheduler_addtask(s, task_type_stars_sort,
+ task_subtype_none, 0, 0, c, NULL);
+ }
+
/* Local tasks only... */
if (c->nodeID == e->nodeID) {
@@ -739,6 +829,10 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
c->hydro.drift = scheduler_addtask(s, task_type_drift_part,
task_subtype_none, 0, 0, c, NULL);
+ /* Add the task finishing the force calculation */
+ c->hydro.end_force = scheduler_addtask(s, task_type_end_hydro_force,
+ task_subtype_none, 0, 0, c, NULL);
+
/* Generate the ghost tasks. */
c->hydro.ghost_in =
scheduler_addtask(s, task_type_ghost_in, task_subtype_none, 0,
@@ -748,57 +842,63 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
/* implicit = */ 1, c, NULL);
engine_add_ghosts(e, c, c->hydro.ghost_in, c->hydro.ghost_out);
-#ifdef EXTRA_HYDRO_LOOP
/* Generate the extra ghost task. */
+#ifdef EXTRA_HYDRO_LOOP
c->hydro.extra_ghost = scheduler_addtask(
s, task_type_extra_ghost, task_subtype_none, 0, 0, c, NULL);
#endif
- }
- } else { /* We are above the super-cell so need to go deeper */
+ /* Stars */
+ if (with_stars) {
+ c->stars.drift = scheduler_addtask(s, task_type_drift_spart,
+ task_subtype_none, 0, 0, c, NULL);
+ scheduler_addunlock(s, c->stars.drift, c->super->kick2);
+ }
- /* Recurse. */
- if (c->split)
- for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL)
- engine_make_hierarchical_tasks_hydro(e, c->progeny[k]);
- }
-}
+ /* Subgrid tasks: cooling */
+ if (with_cooling) {
-/**
- * @brief Generate the stars hierarchical tasks for a hierarchy of cells -
- * i.e. all the O(Npart) tasks -- star version
- *
- * Tasks are only created here. The dependencies will be added later on.
- *
- * Note that there is no need to recurse below the super-cell. Note also
- * that we only add tasks if the relevant particles are present in the cell.
- *
- * @param e The #engine.
- * @param c The #cell.
- */
-void engine_make_hierarchical_tasks_stars(struct engine *e, struct cell *c) {
+ c->hydro.cooling = scheduler_addtask(s, task_type_cooling,
+ task_subtype_none, 0, 0, c, NULL);
- struct scheduler *s = &e->sched;
+ scheduler_addunlock(s, c->hydro.end_force, c->hydro.cooling);
+ scheduler_addunlock(s, c->hydro.cooling, c->super->kick2);
- /* Are we in a super-cell ? */
- if (c->super == c) {
+ } else {
+ scheduler_addunlock(s, c->hydro.end_force, c->super->kick2);
+ }
- /* Add the sort task. */
- c->stars.sorts = scheduler_addtask(s, task_type_stars_sort,
- task_subtype_none, 0, 0, c, NULL);
+ /* Subgrid tasks: star formation */
+ if (with_star_formation) {
- /* Local tasks only... */
- if (c->nodeID == e->nodeID) {
+ c->hydro.star_formation = scheduler_addtask(
+ s, task_type_star_formation, task_subtype_none, 0, 0, c, NULL);
- /* Generate the ghost tasks. */
- c->stars.ghost_in =
- scheduler_addtask(s, task_type_stars_ghost_in, task_subtype_none, 0,
- /* implicit = */ 1, c, NULL);
- c->stars.ghost_out =
- scheduler_addtask(s, task_type_stars_ghost_out, task_subtype_none, 0,
- /* implicit = */ 1, c, NULL);
- engine_add_stars_ghosts(e, c, c->stars.ghost_in, c->stars.ghost_out);
+ scheduler_addunlock(s, c->super->kick2, c->hydro.star_formation);
+ scheduler_addunlock(s, c->hydro.star_formation, c->super->timestep);
+ }
+
+ /* Subgrid tasks: feedback */
+ if (with_feedback) {
+
+ c->stars.stars_in =
+ scheduler_addtask(s, task_type_stars_in, task_subtype_none, 0,
+ /* implicit = */ 1, c, NULL);
+
+ c->stars.stars_out =
+ scheduler_addtask(s, task_type_stars_out, task_subtype_none, 0,
+ /* implicit = */ 1, c, NULL);
+
+ c->stars.ghost = scheduler_addtask(s, task_type_stars_ghost,
+ task_subtype_none, 0, 0, c, NULL);
+
+ scheduler_addunlock(s, c->super->kick2, c->stars.stars_in);
+ scheduler_addunlock(s, c->stars.stars_out, c->super->timestep);
+
+ if (with_star_formation) {
+ scheduler_addunlock(s, c->hydro.star_formation, c->stars.stars_in);
+ }
+ }
}
} else { /* We are above the super-cell so need to go deeper */
@@ -806,7 +906,27 @@ void engine_make_hierarchical_tasks_stars(struct engine *e, struct cell *c) {
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL)
- engine_make_hierarchical_tasks_stars(e, c->progeny[k]);
+ engine_make_hierarchical_tasks_hydro(e, c->progeny[k]);
+ }
+}
+
+void engine_make_hierarchical_tasks_mapper(void *map_data, int num_elements,
+ void *extra_data) {
+
+ struct engine *e = (struct engine *)extra_data;
+ const int with_hydro = (e->policy & engine_policy_hydro);
+ const int with_self_gravity = (e->policy & engine_policy_self_gravity);
+ const int with_ext_gravity = (e->policy & engine_policy_external_gravity);
+
+ for (int ind = 0; ind < num_elements; ind++) {
+ struct cell *c = &((struct cell *)map_data)[ind];
+ /* Make the common tasks (time integration) */
+ engine_make_hierarchical_tasks_common(e, c);
+ /* Add the hydro stuff */
+ if (with_hydro) engine_make_hierarchical_tasks_hydro(e, c);
+ /* And the gravity stuff */
+ if (with_self_gravity || with_ext_gravity)
+ engine_make_hierarchical_tasks_gravity(e, c);
}
}
@@ -973,28 +1093,6 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
}
}
-void engine_make_hierarchical_tasks_mapper(void *map_data, int num_elements,
- void *extra_data) {
- struct engine *e = (struct engine *)extra_data;
- const int is_with_hydro = (e->policy & engine_policy_hydro);
- const int is_with_self_gravity = (e->policy & engine_policy_self_gravity);
- const int is_with_external_gravity =
- (e->policy & engine_policy_external_gravity);
- const int is_with_feedback = (e->policy & engine_policy_feedback);
-
- for (int ind = 0; ind < num_elements; ind++) {
- struct cell *c = &((struct cell *)map_data)[ind];
- /* Make the common tasks (time integration) */
- engine_make_hierarchical_tasks_common(e, c);
- /* Add the hydro stuff */
- if (is_with_hydro) engine_make_hierarchical_tasks_hydro(e, c);
- /* And the gravity stuff */
- if (is_with_self_gravity || is_with_external_gravity)
- engine_make_hierarchical_tasks_gravity(e, c);
- if (is_with_feedback) engine_make_hierarchical_tasks_stars(e, c);
- }
-}
-
/**
* @brief Constructs the top-level tasks for the external gravity.
*
@@ -1056,9 +1154,10 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
/* Link stars sort tasks to all the higher sort task. */
if (t_type == task_type_stars_sort) {
for (struct cell *finger = t->ci->parent; finger != NULL;
- finger = finger->parent)
+ finger = finger->parent) {
if (finger->stars.sorts != NULL)
scheduler_addunlock(sched, t, finger->stars.sorts);
+ }
}
/* Link self tasks to cells. */
@@ -1220,9 +1319,9 @@ void engine_link_gravity_tasks(struct engine *e) {
if (ci_nodeID != nodeID) error("Non-local self task");
#endif
- /* drift -----> gravity --> end_force */
+ /* drift -----> gravity --> end_gravity_force */
scheduler_addunlock(sched, ci->grav.super->grav.drift, t);
- scheduler_addunlock(sched, t, ci->super->end_force);
+ scheduler_addunlock(sched, t, ci->grav.super->grav.end_force);
}
/* Otherwise, pair interaction? */
@@ -1271,7 +1370,7 @@ void engine_link_gravity_tasks(struct engine *e) {
/* drift -----> gravity --> end_force */
scheduler_addunlock(sched, ci->grav.super->grav.drift, t);
- scheduler_addunlock(sched, t, ci->super->end_force);
+ scheduler_addunlock(sched, t, ci->grav.super->grav.end_force);
}
/* Otherwise, sub-pair interaction? */
@@ -1368,22 +1467,6 @@ static inline void engine_make_hydro_loops_dependencies(
}
#endif
-/**
- * @brief Creates the dependency network for the stars tasks of a given cell.
- *
- * @param sched The #scheduler.
- * @param density The star density task to link.
- * @param feedback The star feedback task to link.
- * @param c The cell.
- */
-static inline void engine_make_stars_loops_dependencies(struct scheduler *sched,
- struct task *density,
- struct task *feedback,
- struct cell *c) {
- /* density loop --> ghost --> feedback loop*/
- scheduler_addunlock(sched, density, c->super->stars.ghost_in);
- scheduler_addunlock(sched, c->super->stars.ghost_out, feedback);
-}
/**
* @brief Duplicates the first hydro loop and construct all the
@@ -1403,638 +1486,544 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
const int nodeID = e->nodeID;
const int with_cooling = (e->policy & engine_policy_cooling);
const int with_limiter = (e->policy & engine_policy_limiter);
+ const int with_feedback = (e->policy & engine_policy_feedback);
#ifdef EXTRA_HYDRO_LOOP
struct task *t_gradient = NULL;
#endif
struct task *t_force = NULL;
struct task *t_limiter = NULL;
+ struct task *t_star_density = NULL;
+ struct task *t_star_feedback = NULL;
for (int ind = 0; ind < num_elements; ind++) {
+
struct task *t = &((struct task *)map_data)[ind];
+ const enum task_types t_type = t->type;
+ const enum task_subtypes t_subtype = t->subtype;
+ const long long flags = t->flags;
+ struct cell *ci = t->ci;
+ struct cell *cj = t->cj;
+
+ /* Sort tasks depend on the drift of the cell (gas version). */
+ if (t_type == task_type_sort && ci->nodeID == nodeID) {
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t);
+ }
- /* Sort tasks depend on the drift of the cell. */
- if (t->type == task_type_sort && t->ci->nodeID == engine_rank) {
- scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
+ /* Sort tasks depend on the drift of the cell (stars version). */
+ else if (t_type == task_type_stars_sort && ci->nodeID == nodeID) {
+ scheduler_addunlock(sched, ci->hydro.super->stars.drift, t);
}
/* Self-interaction? */
- else if (t->type == task_type_self && t->subtype == task_subtype_density) {
+ else if (t_type == task_type_self && t_subtype == task_subtype_density) {
/* Make the self-density tasks depend on the drift only. */
- scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t);
-#ifdef EXTRA_HYDRO_LOOP
- /* Start by constructing the task for the second and third hydro loop. */
- t_gradient = scheduler_addtask(sched, task_type_self,
- task_subtype_gradient, 0, 0, t->ci, NULL);
- t_force = scheduler_addtask(sched, task_type_self, task_subtype_force, 0,
- 0, t->ci, NULL);
+ /* Task for the second hydro loop, */
+ t_force = scheduler_addtask(sched, task_type_self, task_subtype_force,
+ flags, 0, ci, NULL);
- /* and the task for the time-step limiter */
- if (with_limiter)
+ /* the task for the time-step limiter */
+ if (with_limiter) {
t_limiter = scheduler_addtask(sched, task_type_self,
- task_subtype_limiter, 0, 0, t->ci, NULL);
+ task_subtype_limiter, flags, 0, ci, NULL);
+ }
+
+ /* The stellar feedback tasks */
+ if (with_feedback) {
+ t_star_density =
+ scheduler_addtask(sched, task_type_self, task_subtype_stars_density,
+ flags, 0, ci, NULL);
+ t_star_feedback =
+ scheduler_addtask(sched, task_type_self,
+ task_subtype_stars_feedback, flags, 0, ci, NULL);
+ }
+
+ /* Link the tasks to the cells */
+ engine_addlink(e, &ci->hydro.force, t_force);
+ if (with_limiter) {
+ engine_addlink(e, &ci->hydro.limiter, t_limiter);
+ }
+ if (with_feedback) {
+ engine_addlink(e, &ci->stars.density, t_star_density);
+ engine_addlink(e, &ci->stars.feedback, t_star_feedback);
+ }
+
+#ifdef EXTRA_HYDRO_LOOP
+
+ /* Same work for the additional hydro loop */
+ t_gradient = scheduler_addtask(sched, task_type_self,
+ task_subtype_gradient, flags, 0, ci, NULL);
/* Add the link between the new loops and the cell */
- engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
- engine_addlink(e, &t->ci->hydro.force, t_force);
- if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
+ engine_addlink(e, &ci->hydro.gradient, t_gradient);
/* Now, build all the dependencies for the hydro */
engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
- t_limiter, t->ci, with_cooling,
+ t_limiter, ci, with_cooling,
with_limiter);
- scheduler_addunlock(sched, t_force, t->ci->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
#else
- /* Start by constructing the task for the second hydro loop */
- t_force = scheduler_addtask(sched, task_type_self, task_subtype_force, 0,
- 0, t->ci, NULL);
+ /* Now, build all the dependencies for the hydro */
+ engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, ci,
+ with_cooling, with_limiter);
+#endif
- /* and the task for the time-step limiter */
- if (with_limiter)
- t_limiter = scheduler_addtask(sched, task_type_self,
- task_subtype_limiter, 0, 0, t->ci, NULL);
-
- /* Add the link between the new loop and the cell */
- engine_addlink(e, &t->ci->hydro.force, t_force);
- if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
+ /* Create the task dependencies */
+ scheduler_addunlock(sched, t_force, ci->hydro.super->hydro.end_force);
+
+ if (with_feedback) {
+
+ scheduler_addunlock(sched, ci->hydro.super->stars.drift,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->stars.stars_in,
+ t_star_density);
+ scheduler_addunlock(sched, t_star_density,
+ ci->hydro.super->stars.ghost);
+ scheduler_addunlock(sched, ci->hydro.super->stars.ghost,
+ t_star_feedback);
+ scheduler_addunlock(sched, t_star_feedback,
+ ci->hydro.super->stars.stars_out);
+ }
- /* Now, build all the dependencies for the hydro */
- engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, t->ci,
- with_cooling, with_limiter);
- scheduler_addunlock(sched, t_force, t->ci->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
-#endif
+ if (with_limiter) {
+ scheduler_addunlock(sched, ci->super->kick2, t_limiter);
+ scheduler_addunlock(sched, t_limiter, ci->super->timestep);
+ scheduler_addunlock(sched, t_limiter, ci->super->timestep_limiter);
+ }
}
/* Otherwise, pair interaction? */
- else if (t->type == task_type_pair && t->subtype == task_subtype_density) {
+ else if (t_type == task_type_pair && t_subtype == task_subtype_density) {
- /* Make all density tasks depend on the drift and the sorts. */
- if (t->ci->nodeID == engine_rank)
- scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
- scheduler_addunlock(sched, t->ci->hydro.super->hydro.sorts, t);
- if (t->ci->hydro.super != t->cj->hydro.super) {
- if (t->cj->nodeID == engine_rank)
- scheduler_addunlock(sched, t->cj->hydro.super->hydro.drift, t);
- scheduler_addunlock(sched, t->cj->hydro.super->hydro.sorts, t);
+ /* Make all density tasks depend on the drift */
+ if (ci->nodeID == nodeID) {
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t);
+ }
+ if ((cj->nodeID == nodeID) && (ci->hydro.super != cj->hydro.super)) {
+ scheduler_addunlock(sched, cj->hydro.super->hydro.drift, t);
}
-#ifdef EXTRA_HYDRO_LOOP
- /* Start by constructing the task for the second and third hydro loop */
- t_gradient = scheduler_addtask(sched, task_type_pair,
- task_subtype_gradient, 0, 0, t->ci, t->cj);
- t_force = scheduler_addtask(sched, task_type_pair, task_subtype_force, 0,
- 0, t->ci, t->cj);
+ /* Make all density tasks depend on the sorts */
+ scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t);
+ if (ci->hydro.super != cj->hydro.super) {
+ scheduler_addunlock(sched, cj->hydro.super->hydro.sorts, t);
+ }
+
+ /* New task for the force */
+ t_force = scheduler_addtask(sched, task_type_pair, task_subtype_force,
+ flags, 0, ci, cj);
/* and the task for the time-step limiter */
- if (with_limiter)
+ if (with_limiter) {
t_limiter = scheduler_addtask(sched, task_type_pair,
- task_subtype_limiter, 0, 0, t->ci, t->cj);
+ task_subtype_limiter, flags, 0, ci, cj);
+ }
+
+ /* The stellar feedback tasks */
+ if (with_feedback) {
+ t_star_density =
+ scheduler_addtask(sched, task_type_pair, task_subtype_stars_density,
+ flags, 0, ci, cj);
+ t_star_feedback =
+ scheduler_addtask(sched, task_type_pair,
+ task_subtype_stars_feedback, flags, 0, ci, cj);
+ }
+
+ engine_addlink(e, &ci->hydro.force, t_force);
+ engine_addlink(e, &cj->hydro.force, t_force);
+ if (with_limiter) {
+ engine_addlink(e, &ci->hydro.limiter, t_limiter);
+ engine_addlink(e, &cj->hydro.limiter, t_limiter);
+ }
+ if (with_feedback) {
+ engine_addlink(e, &ci->stars.density, t_star_density);
+ engine_addlink(e, &cj->stars.density, t_star_density);
+ engine_addlink(e, &ci->stars.feedback, t_star_feedback);
+ engine_addlink(e, &cj->stars.feedback, t_star_feedback);
+ }
+
+#ifdef EXTRA_HYDRO_LOOP
+
+ /* Start by constructing the task for the second and third hydro loop */
+ t_gradient = scheduler_addtask(sched, task_type_pair,
+ task_subtype_gradient, flags, 0, ci, cj);
/* Add the link between the new loop and both cells */
- engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
- engine_addlink(e, &t->cj->hydro.gradient, t_gradient);
- engine_addlink(e, &t->ci->hydro.force, t_force);
- engine_addlink(e, &t->cj->hydro.force, t_force);
- if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
- if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
+ engine_addlink(e, &ci->hydro.gradient, t_gradient);
+ engine_addlink(e, &cj->hydro.gradient, t_gradient);
/* Now, build all the dependencies for the hydro for the cells */
/* that are local and are not descendant of the same super_hydro-cells */
- if (t->ci->nodeID == nodeID) {
+ if (ci->nodeID == nodeID) {
engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
- t_limiter, t->ci, with_cooling,
+ t_limiter, ci, with_cooling,
+ with_limiter);
+ }
+ if ((cj->nodeID == nodeID) && (ci->hydro.super != cj->hydro.super)) {
+ engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
+ t_limiter, cj, with_cooling,
with_limiter);
- scheduler_addunlock(sched, t_force, t->ci->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
- }
- if (t->cj->nodeID == nodeID) {
- if (t->ci->hydro.super != t->cj->hydro.super) {
- engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
- t_limiter, t->cj, with_cooling,
- with_limiter);
- }
-
- if (t->ci->super != t->cj->super) {
- scheduler_addunlock(sched, t_force, t->cj->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter,
- t->cj->super->timestep_limiter);
- }
}
-
#else
- /* Start by constructing the task for the second hydro loop */
- t_force = scheduler_addtask(sched, task_type_pair, task_subtype_force, 0,
- 0, t->ci, t->cj);
+ /* Now, build all the dependencies for the hydro for the cells */
+ /* that are local and are not descendant of the same super_hydro-cells */
+ if (ci->nodeID == nodeID) {
+ engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, ci,
+ with_cooling, with_limiter);
+ }
+ if ((cj->nodeID == nodeID) && (ci->hydro.super != cj->hydro.super)) {
+ engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, cj,
+ with_cooling, with_limiter);
+ }
+#endif
- /* and the task for the time-step limiter */
- if (with_limiter)
- t_limiter = scheduler_addtask(sched, task_type_pair,
- task_subtype_limiter, 0, 0, t->ci, t->cj);
+ if (with_feedback) {
+ scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
+ t_star_density);
- /* Add the link between the new loop and both cells */
- engine_addlink(e, &t->ci->hydro.force, t_force);
- engine_addlink(e, &t->cj->hydro.force, t_force);
- if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
- if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
+ if (ci->hydro.super != cj->hydro.super) {
+ scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
+ t_star_density);
+ }
+ }
- /* Now, build all the dependencies for the hydro for the cells */
- /* that are local and are not descendant of the same super_hydro-cells */
- if (t->ci->nodeID == nodeID) {
- engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter,
- t->ci, with_cooling, with_limiter);
- scheduler_addunlock(sched, t_force, t->ci->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
- }
- if (t->cj->nodeID == nodeID) {
- if (t->ci->hydro.super != t->cj->hydro.super) {
- engine_make_hydro_loops_dependencies(
- sched, t, t_force, t_limiter, t->cj, with_cooling, with_limiter);
+ if (ci->nodeID == nodeID) {
+ scheduler_addunlock(sched, t_force, ci->hydro.super->hydro.end_force);
+
+ if (with_feedback) {
+
+ scheduler_addunlock(sched, ci->hydro.super->stars.drift,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->stars.stars_in,
+ t_star_density);
+ scheduler_addunlock(sched, t_star_density,
+ ci->hydro.super->stars.ghost);
+ scheduler_addunlock(sched, ci->hydro.super->stars.ghost,
+ t_star_feedback);
+ scheduler_addunlock(sched, t_star_feedback,
+ ci->hydro.super->stars.stars_out);
}
- if (t->ci->super != t->cj->super) {
- scheduler_addunlock(sched, t_force, t->cj->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter,
- t->cj->super->timestep_limiter);
+ if (with_limiter) {
+ scheduler_addunlock(sched, ci->super->kick2, t_limiter);
+ scheduler_addunlock(sched, t_limiter, ci->super->timestep);
+ scheduler_addunlock(sched, t_limiter, ci->super->timestep_limiter);
+ }
+ } else /*(ci->nodeID != nodeID) */ {
+ if (with_feedback) {
+ scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
+ t_star_feedback);
}
}
-#endif
+ if (cj->nodeID == nodeID) {
- }
+ if (ci->hydro.super != cj->hydro.super) {
- /* Otherwise, sub-self interaction? */
- else if (t->type == task_type_sub_self &&
- t->subtype == task_subtype_density) {
+ scheduler_addunlock(sched, t_force, cj->hydro.super->hydro.end_force);
- /* Make all density tasks depend on the drift and sorts. */
- scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
- scheduler_addunlock(sched, t->ci->hydro.super->hydro.sorts, t);
-
-#ifdef EXTRA_HYDRO_LOOP
+ if (with_feedback) {
- /* Start by constructing the task for the second and third hydro loop */
- t_gradient =
- scheduler_addtask(sched, task_type_sub_self, task_subtype_gradient,
- t->flags, 0, t->ci, NULL);
- t_force = scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
- t->flags, 0, t->ci, NULL);
+ scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
+ t_star_density);
+ scheduler_addunlock(sched, cj->hydro.super->stars.drift,
+ t_star_density);
+ scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
+ t_star_density);
+ scheduler_addunlock(sched, cj->hydro.super->stars.stars_in,
+ t_star_density);
+ scheduler_addunlock(sched, t_star_density,
+ cj->hydro.super->stars.ghost);
+ scheduler_addunlock(sched, cj->hydro.super->stars.ghost,
+ t_star_feedback);
+ scheduler_addunlock(sched, t_star_feedback,
+ cj->hydro.super->stars.stars_out);
+ }
- /* and the task for the time-step limiter */
- if (with_limiter)
- t_limiter =
- scheduler_addtask(sched, task_type_sub_self, task_subtype_limiter,
- t->flags, 0, t->ci, NULL);
+ if (with_limiter) {
+ scheduler_addunlock(sched, cj->super->kick2, t_limiter);
+ scheduler_addunlock(sched, t_limiter, cj->super->timestep);
+ scheduler_addunlock(sched, t_limiter, cj->super->timestep_limiter);
+ }
+ }
+ } else /*(cj->nodeID != nodeID) */ {
+ if (with_feedback) {
+ scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
+ t_star_feedback);
+ }
+ }
+ }
- /* Add the link between the new loop and the cell */
- engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
- engine_addlink(e, &t->ci->hydro.force, t_force);
- if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
+ /* Otherwise, sub-self interaction? */
+ else if (t_type == task_type_sub_self &&
+ t_subtype == task_subtype_density) {
- /* Now, build all the dependencies for the hydro for the cells */
- /* that are local and are not descendant of the same super_hydro-cells */
- if (t->ci->nodeID == nodeID) {
- engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
- t_limiter, t->ci, with_cooling,
- with_limiter);
- scheduler_addunlock(sched, t_force, t->ci->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
- }
+ /* Make all density tasks depend on the drift and sorts. */
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t);
-#else
/* Start by constructing the task for the second hydro loop */
t_force = scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
- t->flags, 0, t->ci, NULL);
+ flags, 0, ci, NULL);
/* and the task for the time-step limiter */
- if (with_limiter)
- t_limiter =
- scheduler_addtask(sched, task_type_sub_self, task_subtype_limiter,
- t->flags, 0, t->ci, NULL);
-
- /* Add the link between the new loop and the cell */
- engine_addlink(e, &t->ci->hydro.force, t_force);
- if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
-
- /* Now, build all the dependencies for the hydro for the cells */
- /* that are local and are not descendant of the same super_hydro-cells */
- if (t->ci->nodeID == nodeID) {
- engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter,
- t->ci, with_cooling, with_limiter);
- scheduler_addunlock(sched, t_force, t->ci->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
+ if (with_limiter) {
+ t_limiter = scheduler_addtask(sched, task_type_sub_self,
+ task_subtype_limiter, flags, 0, ci, NULL);
}
-#endif
- }
- /* Otherwise, sub-pair interaction? */
- else if (t->type == task_type_sub_pair &&
- t->subtype == task_subtype_density) {
+ /* The stellar feedback tasks */
+ if (with_feedback) {
+ t_star_density =
+ scheduler_addtask(sched, task_type_sub_self,
+ task_subtype_stars_density, flags, 0, ci, NULL);
+ t_star_feedback =
+ scheduler_addtask(sched, task_type_sub_self,
+ task_subtype_stars_feedback, flags, 0, ci, NULL);
+ }
- /* Make all density tasks depend on the drift. */
- if (t->ci->nodeID == engine_rank)
- scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
- scheduler_addunlock(sched, t->ci->hydro.super->hydro.sorts, t);
- if (t->ci->hydro.super != t->cj->hydro.super) {
- if (t->cj->nodeID == engine_rank)
- scheduler_addunlock(sched, t->cj->hydro.super->hydro.drift, t);
- scheduler_addunlock(sched, t->cj->hydro.super->hydro.sorts, t);
+ /* Add the link between the new loop and the cell */
+ engine_addlink(e, &ci->hydro.force, t_force);
+ if (with_limiter) {
+ engine_addlink(e, &ci->hydro.limiter, t_limiter);
+ }
+ if (with_feedback) {
+ engine_addlink(e, &ci->stars.density, t_star_density);
+ engine_addlink(e, &ci->stars.feedback, t_star_feedback);
}
#ifdef EXTRA_HYDRO_LOOP
/* Start by constructing the task for the second and third hydro loop */
- t_gradient =
- scheduler_addtask(sched, task_type_sub_pair, task_subtype_gradient,
- t->flags, 0, t->ci, t->cj);
- t_force = scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
- t->flags, 0, t->ci, t->cj);
+ t_gradient = scheduler_addtask(sched, task_type_sub_self,
+ task_subtype_gradient, flags, 0, ci, NULL);
- /* and the task for the time-step limiter */
- if (with_limiter)
- t_limiter =
- scheduler_addtask(sched, task_type_sub_pair, task_subtype_limiter,
- t->flags, 0, t->ci, t->cj);
-
- /* Add the link between the new loop and both cells */
- engine_addlink(e, &t->ci->hydro.gradient, t_gradient);
- engine_addlink(e, &t->cj->hydro.gradient, t_gradient);
- engine_addlink(e, &t->ci->hydro.force, t_force);
- engine_addlink(e, &t->cj->hydro.force, t_force);
- if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
- if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
+ /* Add the link between the new loop and the cell */
+ engine_addlink(e, &ci->hydro.gradient, t_gradient);
/* Now, build all the dependencies for the hydro for the cells */
/* that are local and are not descendant of the same super_hydro-cells */
- if (t->ci->nodeID == nodeID) {
- engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
- t_limiter, t->ci, with_cooling,
- with_limiter);
- scheduler_addunlock(sched, t_force, t->ci->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
- }
- if (t->cj->nodeID == nodeID) {
- if (t->ci->hydro.super != t->cj->hydro.super) {
- engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
- t_limiter, t->cj, with_cooling,
- with_limiter);
- }
-
- if (t->ci->super != t->cj->super) {
- scheduler_addunlock(sched, t_force, t->cj->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter,
- t->cj->super->timestep_limiter);
- }
- }
-
+ engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
+ t_limiter, ci, with_cooling,
+ with_limiter);
#else
- /* Start by constructing the task for the second hydro loop */
- t_force = scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
- t->flags, 0, t->ci, t->cj);
-
- /* and the task for the time-step limiter */
- if (with_limiter)
- t_limiter =
- scheduler_addtask(sched, task_type_sub_pair, task_subtype_limiter,
- t->flags, 0, t->ci, t->cj);
-
- /* Add the link between the new loop and both cells */
- engine_addlink(e, &t->ci->hydro.force, t_force);
- engine_addlink(e, &t->cj->hydro.force, t_force);
- if (with_limiter) engine_addlink(e, &t->ci->hydro.limiter, t_limiter);
- if (with_limiter) engine_addlink(e, &t->cj->hydro.limiter, t_limiter);
/* Now, build all the dependencies for the hydro for the cells */
/* that are local and are not descendant of the same super_hydro-cells */
- if (t->ci->nodeID == nodeID) {
- engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter,
- t->ci, with_cooling, with_limiter);
-
- scheduler_addunlock(sched, t_force, t->ci->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->ci->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->ci->super->timestep_limiter);
- }
- if (t->cj->nodeID == nodeID) {
- if (t->ci->hydro.super != t->cj->hydro.super) {
- engine_make_hydro_loops_dependencies(
- sched, t, t_force, t_limiter, t->cj, with_cooling, with_limiter);
- }
-
- if (t->ci->super != t->cj->super) {
- scheduler_addunlock(sched, t_force, t->cj->super->end_force);
- if (with_limiter)
- scheduler_addunlock(sched, t->cj->super->kick2, t_limiter);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter, t->cj->super->timestep);
- if (with_limiter)
- scheduler_addunlock(sched, t_limiter,
- t->cj->super->timestep_limiter);
- }
- }
+ engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, ci,
+ with_cooling, with_limiter);
#endif
- }
- }
-}
-
-/**
- * @brief Duplicates the first stars loop and construct all the
- * dependencies for the stars part
- *
- * This is done by looping over all the previously constructed tasks
- * and adding another task involving the same cells but this time
- * corresponding to the second stars loop over neighbours.
- * With all the relevant tasks for a given cell available, we construct
- * all the dependencies for that cell.
- */
-void engine_make_extra_starsloop_tasks_mapper(void *map_data, int num_elements,
- void *extra_data) {
- struct engine *e = (struct engine *)extra_data;
- struct scheduler *sched = &e->sched;
- const int nodeID = e->nodeID;
-
- for (int ind = 0; ind < num_elements; ind++) {
- struct task *t = &((struct task *)map_data)[ind];
-
- /* Sort tasks depend on the drift and gravity drift of the cell. */
- if (t->type == task_type_stars_sort && t->ci->nodeID == engine_rank) {
- scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
- scheduler_addunlock(sched, t->ci->super->grav.drift, t);
- }
-
- /* Self-interaction? */
- else if (t->type == task_type_self &&
- t->subtype == task_subtype_stars_density) {
-
- /* Make the self-density tasks depend on the drift and gravity drift. */
- scheduler_addunlock(sched, t->ci->hydro.super->hydro.drift, t);
- scheduler_addunlock(sched, t->ci->super->grav.drift, t);
-
- /* Start by constructing the task for the second stars loop */
- struct task *t2 =
- scheduler_addtask(sched, task_type_self, task_subtype_stars_feedback,
- 0, 0, t->ci, NULL);
-
- /* Add the link between the new loop and the cell */
- engine_addlink(e, &t->ci->stars.feedback, t2);
+ /* Create the task dependencies */
+ scheduler_addunlock(sched, t_force, ci->hydro.super->hydro.end_force);
+
+ if (with_feedback) {
+
+ scheduler_addunlock(sched, ci->hydro.super->stars.drift,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->stars.stars_in,
+ t_star_density);
+ scheduler_addunlock(sched, t_star_density,
+ ci->hydro.super->stars.ghost);
+ scheduler_addunlock(sched, ci->hydro.super->stars.ghost,
+ t_star_feedback);
+ scheduler_addunlock(sched, t_star_feedback,
+ ci->hydro.super->stars.stars_out);
+ }
- /* Now, build all the dependencies for the stars */
- engine_make_stars_loops_dependencies(sched, t, t2, t->ci);
+ if (with_limiter) {
+ scheduler_addunlock(sched, ci->super->kick2, t_limiter);
+ scheduler_addunlock(sched, t_limiter, ci->super->timestep);
+ scheduler_addunlock(sched, t_limiter, ci->super->timestep_limiter);
+ }
- /* end_force depends on feedback tasks */
- scheduler_addunlock(sched, t2, t->ci->super->end_force);
}
- /* Otherwise, pair interaction? */
- else if (t->type == task_type_pair &&
- t->subtype == task_subtype_stars_density) {
-
- /* Make all stars density tasks depend on the hydro drift and sorts,
- * gravity drift and star sorts. */
- if (t->ci->nodeID == engine_rank)
- scheduler_addunlock(sched, t->ci->super->hydro.drift, t);
- scheduler_addunlock(sched, t->ci->super->hydro.sorts, t);
- if (t->cj->nodeID == engine_rank)
- scheduler_addunlock(sched, t->cj->super->grav.drift, t);
- scheduler_addunlock(sched, t->ci->super->stars.sorts, t);
-
- if (t->ci->super != t->cj->super) {
- if (t->cj->nodeID == engine_rank)
- scheduler_addunlock(sched, t->cj->super->hydro.drift, t);
- scheduler_addunlock(sched, t->cj->super->hydro.sorts, t);
- if (t->ci->nodeID == engine_rank)
- scheduler_addunlock(sched, t->ci->super->grav.drift, t);
- scheduler_addunlock(sched, t->cj->super->stars.sorts, t);
- }
-
- /* Start by constructing the task for the second stars loop */
- struct task *t2 =
- scheduler_addtask(sched, task_type_pair, task_subtype_stars_feedback,
- 0, 0, t->ci, t->cj);
-
- /* Add the link between the new loop and both cells */
- engine_addlink(e, &t->ci->stars.feedback, t2);
- engine_addlink(e, &t->cj->stars.feedback, t2);
+ /* Otherwise, sub-pair interaction? */
+ else if (t_type == task_type_sub_pair &&
+ t_subtype == task_subtype_density) {
- /* Now, build all the dependencies for the stars for the cells */
- if (t->ci->nodeID == nodeID) {
- engine_make_stars_loops_dependencies(sched, t, t2, t->ci);
- scheduler_addunlock(sched, t2, t->ci->super->end_force);
+ /* Make all density tasks depend on the drift */
+ if (ci->nodeID == nodeID) {
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t);
}
- if (t->cj->nodeID == nodeID) {
- if (t->ci->super != t->cj->super)
- engine_make_stars_loops_dependencies(sched, t, t2, t->cj);
- if (t->ci->super != t->cj->super)
- scheduler_addunlock(sched, t2, t->cj->super->end_force);
+ if ((cj->nodeID == nodeID) && (ci->hydro.super != cj->hydro.super)) {
+ scheduler_addunlock(sched, cj->hydro.super->hydro.drift, t);
}
- }
- /* Otherwise, sub-self interaction? */
- else if (t->type == task_type_sub_self &&
- t->subtype == task_subtype_stars_density) {
+ /* Make all density tasks depend on the sorts */
+ scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t);
+ if (ci->hydro.super != cj->hydro.super) {
+ scheduler_addunlock(sched, cj->hydro.super->hydro.sorts, t);
+ }
- /* Make all stars density tasks depend on the hydro drift and sorts,
- * gravity drift and star sorts. */
- scheduler_addunlock(sched, t->ci->super->hydro.drift, t);
- scheduler_addunlock(sched, t->ci->super->hydro.sorts, t);
- scheduler_addunlock(sched, t->ci->super->grav.drift, t);
- scheduler_addunlock(sched, t->ci->super->stars.sorts, t);
+ /* New task for the force */
+ t_force = scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
+ flags, 0, ci, cj);
- /* Start by constructing the task for the second stars loop */
- struct task *t2 = scheduler_addtask(sched, task_type_sub_self,
- task_subtype_stars_feedback, t->flags,
- 0, t->ci, t->cj);
+ /* and the task for the time-step limiter */
+ if (with_limiter) {
+ t_limiter = scheduler_addtask(sched, task_type_sub_pair,
+ task_subtype_limiter, flags, 0, ci, cj);
+ }
- /* Add the link between the new loop and the cell */
- engine_addlink(e, &t->ci->stars.feedback, t2);
+ /* The stellar feedback tasks */
+ if (with_feedback) {
+ t_star_density =
+ scheduler_addtask(sched, task_type_sub_pair,
+ task_subtype_stars_density, flags, 0, ci, cj);
+ t_star_feedback =
+ scheduler_addtask(sched, task_type_sub_pair,
+ task_subtype_stars_feedback, flags, 0, ci, cj);
+ }
- /* Now, build all the dependencies for the stars for the cells */
- if (t->ci->nodeID == nodeID) {
- engine_make_stars_loops_dependencies(sched, t, t2, t->ci);
- scheduler_addunlock(sched, t2, t->ci->super->end_force);
+ engine_addlink(e, &ci->hydro.force, t_force);
+ engine_addlink(e, &cj->hydro.force, t_force);
+ if (with_limiter) {
+ engine_addlink(e, &ci->hydro.limiter, t_limiter);
+ engine_addlink(e, &cj->hydro.limiter, t_limiter);
+ }
+ if (with_feedback) {
+ engine_addlink(e, &ci->stars.density, t_star_density);
+ engine_addlink(e, &cj->stars.density, t_star_density);
+ engine_addlink(e, &ci->stars.feedback, t_star_feedback);
+ engine_addlink(e, &cj->stars.feedback, t_star_feedback);
}
- }
- /* Otherwise, sub-pair interaction? */
- else if (t->type == task_type_sub_pair &&
- t->subtype == task_subtype_stars_density) {
-
- /* Make all stars density tasks depend on the hydro drift and sorts,
- * gravity drift and star sorts. */
- if (t->cj->nodeID == engine_rank)
- scheduler_addunlock(sched, t->cj->super->hydro.drift, t);
- scheduler_addunlock(sched, t->cj->super->hydro.sorts, t);
- if (t->cj->nodeID == engine_rank)
- scheduler_addunlock(sched, t->cj->super->grav.drift, t);
- scheduler_addunlock(sched, t->ci->super->stars.sorts, t);
-
- if (t->ci->super != t->cj->super) {
- if (t->ci->nodeID == engine_rank)
- scheduler_addunlock(sched, t->ci->super->hydro.drift, t);
- scheduler_addunlock(sched, t->ci->super->hydro.sorts, t);
- if (t->ci->nodeID == engine_rank)
- scheduler_addunlock(sched, t->ci->super->grav.drift, t);
- scheduler_addunlock(sched, t->cj->super->stars.sorts, t);
- }
-
- /* Start by constructing the task for the second stars loop */
- struct task *t2 = scheduler_addtask(sched, task_type_sub_pair,
- task_subtype_stars_feedback, t->flags,
- 0, t->ci, t->cj);
+#ifdef EXTRA_HYDRO_LOOP
+
+ /* Start by constructing the task for the second and third hydro loop */
+ t_gradient = scheduler_addtask(sched, task_type_sub_pair,
+ task_subtype_gradient, flags, 0, ci, cj);
/* Add the link between the new loop and both cells */
- engine_addlink(e, &t->ci->stars.feedback, t2);
- engine_addlink(e, &t->cj->stars.feedback, t2);
+ engine_addlink(e, &ci->hydro.gradient, t_gradient);
+ engine_addlink(e, &cj->hydro.gradient, t_gradient);
- /* Now, build all the dependencies for the stars for the cells */
- if (t->ci->nodeID == nodeID) {
- engine_make_stars_loops_dependencies(sched, t, t2, t->ci);
- scheduler_addunlock(sched, t2, t->ci->super->end_force);
+ /* Now, build all the dependencies for the hydro for the cells */
+ /* that are local and are not descendant of the same super_hydro-cells */
+ if (ci->nodeID == nodeID) {
+ engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
+ t_limiter, ci, with_cooling,
+ with_limiter);
}
- if (t->cj->nodeID == nodeID) {
- if (t->ci->super != t->cj->super)
- engine_make_stars_loops_dependencies(sched, t, t2, t->cj);
- if (t->ci->super != t->cj->super)
- scheduler_addunlock(sched, t2, t->cj->super->end_force);
+ if ((cj->nodeID == nodeID) && (ci->hydro.super != cj->hydro.super)) {
+ engine_make_hydro_loops_dependencies(sched, t, t_gradient, t_force,
+ t_limiter, cj, with_cooling,
+ with_limiter);
}
- }
- }
-}
-
-/**
- * @brief Constructs the top-level pair tasks for the star loop over
- * neighbours
- *
- * Here we construct all the tasks for all possible neighbouring non-empty
- * local cells in the hierarchy. No dependencies are being added thus far.
- * Additional loop over neighbours can later be added by simply duplicating
- * all the tasks created by this function.
- *
- * @param map_data Offset of first two indices disguised as a pointer.
- * @param num_elements Number of cells to traverse.
- * @param extra_data The #engine.
- */
-void engine_make_starsloop_tasks_mapper(void *map_data, int num_elements,
- void *extra_data) {
+#else
- /* Extract the engine pointer. */
- struct engine *e = (struct engine *)extra_data;
- const int periodic = e->s->periodic;
+ /* Now, build all the dependencies for the hydro for the cells */
+ /* that are local and are not descendant of the same super_hydro-cells */
+ if (ci->nodeID == nodeID) {
+ engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, ci,
+ with_cooling, with_limiter);
+ }
+ if ((cj->nodeID == nodeID) && (ci->hydro.super != cj->hydro.super)) {
+ engine_make_hydro_loops_dependencies(sched, t, t_force, t_limiter, cj,
+ with_cooling, with_limiter);
+ }
+#endif
- struct space *s = e->s;
- struct scheduler *sched = &e->sched;
- const int nodeID = e->nodeID;
- const int *cdim = s->cdim;
- struct cell *cells = s->cells_top;
+ if (with_feedback) {
+ scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
+ t_star_density);
+ if (ci->hydro.super != cj->hydro.super) {
+ scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
+ t_star_density);
+ }
+ }
- /* Loop through the elements, which are just byte offsets from NULL. */
- for (int ind = 0; ind < num_elements; ind++) {
+ if (ci->nodeID == nodeID) {
+ scheduler_addunlock(sched, t_force, ci->hydro.super->hydro.end_force);
+
+ if (with_feedback) {
+
+ scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->stars.drift,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+ t_star_density);
+ scheduler_addunlock(sched, ci->hydro.super->stars.stars_in,
+ t_star_density);
+ scheduler_addunlock(sched, t_star_density,
+ ci->hydro.super->stars.ghost);
+ scheduler_addunlock(sched, ci->hydro.super->stars.ghost,
+ t_star_feedback);
+ scheduler_addunlock(sched, t_star_feedback,
+ ci->hydro.super->stars.stars_out);
+ }
- /* Get the cell index. */
- const int cid = (size_t)(map_data) + ind;
- const int i = cid / (cdim[1] * cdim[2]);
- const int j = (cid / cdim[2]) % cdim[1];
- const int k = cid % cdim[2];
+ if (with_limiter) {
+ scheduler_addunlock(sched, ci->super->kick2, t_limiter);
+ scheduler_addunlock(sched, t_limiter, ci->super->timestep);
+ scheduler_addunlock(sched, t_limiter, ci->super->timestep_limiter);
+ }
+ } else /* ci->nodeID != nodeID */ {
- /* Get the cell */
- struct cell *ci = &cells[cid];
+ if (with_feedback) {
+ /* message("%p/%p",ci->hydro.super->stars.sorts, t_star_feedback); */
+ scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
+ t_star_feedback);
+ }
+ }
- /* Skip cells without particles */
- if (ci->stars.count == 0 && ci->hydro.count == 0) continue;
+ if (cj->nodeID == nodeID) {
- /* If the cells is local build a self-interaction */
- if (ci->nodeID == nodeID) {
- scheduler_addtask(sched, task_type_self, task_subtype_stars_density, 0, 0,
- ci, NULL);
- }
+ if (ci->hydro.super != cj->hydro.super) {
- /* Now loop over all the neighbours of this cell */
- for (int ii = -1; ii < 2; ii++) {
- int iii = i + ii;
- if (!periodic && (iii < 0 || iii >= cdim[0])) continue;
- iii = (iii + cdim[0]) % cdim[0];
- for (int jj = -1; jj < 2; jj++) {
- int jjj = j + jj;
- if (!periodic && (jjj < 0 || jjj >= cdim[1])) continue;
- jjj = (jjj + cdim[1]) % cdim[1];
- for (int kk = -1; kk < 2; kk++) {
- int kkk = k + kk;
- if (!periodic && (kkk < 0 || kkk >= cdim[2])) continue;
- kkk = (kkk + cdim[2]) % cdim[2];
+ scheduler_addunlock(sched, t_force, cj->hydro.super->hydro.end_force);
- /* Get the neighbouring cell */
- const int cjd = cell_getid(cdim, iii, jjj, kkk);
- struct cell *cj = &cells[cjd];
+ if (with_feedback) {
- /* Is that neighbour local and does it have particles ? */
- if (cid >= cjd || (cj->stars.count == 0 && cj->hydro.count == 0) ||
- (ci->nodeID != nodeID && cj->nodeID != nodeID))
- continue;
+ scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
+ t_star_density);
+ scheduler_addunlock(sched, cj->hydro.super->stars.drift,
+ t_star_density);
+ scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
+ t_star_density);
+ scheduler_addunlock(sched, cj->hydro.super->stars.stars_in,
+ t_star_density);
+ scheduler_addunlock(sched, t_star_density,
+ cj->hydro.super->stars.ghost);
+ scheduler_addunlock(sched, cj->hydro.super->stars.ghost,
+ t_star_feedback);
+ scheduler_addunlock(sched, t_star_feedback,
+ cj->hydro.super->stars.stars_out);
+ }
- /* Construct the pair task */
- const int sid = sortlistID[(kk + 1) + 3 * ((jj + 1) + 3 * (ii + 1))];
- scheduler_addtask(sched, task_type_pair, task_subtype_stars_density,
- sid, 0, ci, cj);
+ if (with_limiter) {
+ scheduler_addunlock(sched, cj->super->kick2, t_limiter);
+ scheduler_addunlock(sched, t_limiter, cj->super->timestep);
+ scheduler_addunlock(sched, t_limiter, cj->super->timestep_limiter);
+ }
+ }
+ } else /* cj->nodeID != nodeID */ {
+ if (with_feedback) {
+ scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
+ t_star_feedback);
}
}
}
}
}
-
/**
* @brief Constructs the top-level pair tasks for the first hydro loop over
* neighbours
@@ -2054,6 +2043,7 @@ void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* Extract the engine pointer. */
struct engine *e = (struct engine *)extra_data;
const int periodic = e->s->periodic;
+ const int with_feedback = (e->policy & engine_policy_feedback);
struct space *s = e->s;
struct scheduler *sched = &e->sched;
@@ -2075,8 +2065,9 @@ void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* Get the cell */
struct cell *ci = &cells[cid];
- /* Skip cells without hydro particles */
- if (ci->hydro.count == 0) continue;
+ /* Skip cells without hydro or star particles */
+ if ((ci->hydro.count == 0) && (with_feedback && ci->stars.count == 0))
+ continue;
/* If the cell is local build a self-interaction */
if (ci->nodeID == nodeID) {
@@ -2102,8 +2093,10 @@ void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
const int cjd = cell_getid(cdim, iii, jjj, kkk);
struct cell *cj = &cells[cjd];
- /* Is that neighbour local and does it have particles ? */
- if (cid >= cjd || cj->hydro.count == 0 ||
+ /* Is that neighbour local and does it have gas or star particles ? */
+ if ((cid >= cjd) ||
+ ((cj->hydro.count == 0) &&
+ (with_feedback && cj->stars.count == 0)) ||
(ci->nodeID != nodeID && cj->nodeID != nodeID))
continue;
@@ -2168,6 +2161,7 @@ struct cell_type_pair {
void engine_addtasks_send_mapper(void *map_data, int num_elements,
void *extra_data) {
+
struct engine *e = (struct engine *)extra_data;
const int with_limiter = (e->policy & engine_policy_limiter);
struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
@@ -2186,6 +2180,11 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
engine_addtasks_send_hydro(e, ci, cj, /*t_xv=*/NULL,
/*t_rho=*/NULL, /*t_gradient=*/NULL);
+ /* Add the send tasks for the cells in the proxy that have a stars
+ * connection. */
+ if ((e->policy & engine_policy_feedback) && (type & proxy_cell_type_hydro))
+ engine_addtasks_send_stars(e, ci, cj, /*t_feedback=*/NULL);
+
/* Add the send tasks for the cells in the proxy that have a gravity
* connection. */
if ((e->policy & engine_policy_self_gravity) &&
@@ -2196,6 +2195,7 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
void engine_addtasks_recv_mapper(void *map_data, int num_elements,
void *extra_data) {
+
struct engine *e = (struct engine *)extra_data;
const int with_limiter = (e->policy & engine_policy_limiter);
struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
@@ -2212,6 +2212,11 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
if ((e->policy & engine_policy_hydro) && (type & proxy_cell_type_hydro))
engine_addtasks_recv_hydro(e, ci, NULL, NULL, NULL);
+ /* Add the recv tasks for the cells in the proxy that have a stars
+ * connection. */
+ if ((e->policy & engine_policy_feedback) && (type & proxy_cell_type_hydro))
+ engine_addtasks_recv_stars(e, ci, NULL);
+
/* Add the recv tasks for the cells in the proxy that have a gravity
* connection. */
if ((e->policy & engine_policy_self_gravity) &&
@@ -2249,18 +2254,6 @@ void engine_maketasks(struct engine *e) {
tic2 = getticks();
- /* Construct the stars hydro loop over neighbours */
- if (e->policy & engine_policy_feedback) {
- threadpool_map(&e->threadpool, engine_make_starsloop_tasks_mapper, NULL,
- s->nr_cells, 1, 0, e);
- }
-
- if (e->verbose)
- message("Making stellar feedback tasks took %.3f %s.",
- clocks_from_ticks(getticks() - tic2), clocks_getunit());
-
- tic2 = getticks();
-
/* Add the self gravity tasks. */
if (e->policy & engine_policy_self_gravity) {
threadpool_map(&e->threadpool, engine_make_self_gravity_tasks_mapper, NULL,
@@ -2357,19 +2350,6 @@ void engine_maketasks(struct engine *e) {
tic2 = getticks();
- /* Run through the tasks and make stars feedback tasks for each stars density
- task. Each stars feedback task depends on the stars ghosts and unlocks the
- kick task of its super-cell. */
- if (e->policy & engine_policy_stars)
- threadpool_map(&e->threadpool, engine_make_extra_starsloop_tasks_mapper,
- sched->tasks, sched->nr_tasks, sizeof(struct task), 0, e);
-
- if (e->verbose)
- message("Making extra starsloop tasks took %.3f %s.",
- clocks_from_ticks(getticks() - tic2), clocks_getunit());
-
- tic2 = getticks();
-
/* Add the dependencies for the gravity stuff */
if (e->policy & (engine_policy_self_gravity | engine_policy_external_gravity))
engine_link_gravity_tasks(e);
@@ -2381,9 +2361,6 @@ void engine_maketasks(struct engine *e) {
tic2 = getticks();
#ifdef WITH_MPI
- if (e->policy & engine_policy_feedback)
- error("Cannot run stellar feedback with MPI (yet).");
-
/* Add the communication tasks if MPI is being used. */
if (e->policy & engine_policy_mpi) {
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index 3a26dbb2f47f9503aa0b93fa28d679f5eebaeede..c02eb5d2bd272111808701269faed07cef505449 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -84,7 +84,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Local pointer. */
struct cell *ci = t->ci;
- if (ci->nodeID != engine_rank) error("Non-local self task found");
+ if (ci->nodeID != nodeID) error("Non-local self task found");
/* Activate the hydro drift */
if (t_type == task_type_self && t_subtype == task_subtype_density) {
@@ -124,7 +124,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
}
-#ifdef EXTRA_HYDRO_LOOP
else if (t_type == task_type_self && t_subtype == task_subtype_gradient) {
if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
}
@@ -133,7 +132,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
t_subtype == task_subtype_gradient) {
if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
}
-#endif
/* Activate the star density */
else if (t_type == task_type_self &&
@@ -154,7 +152,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
}
- /* Activate the star feedback */
else if (t_type == task_type_self &&
t_subtype == task_subtype_stars_feedback) {
if (cell_is_active_stars(ci, e)) {
@@ -162,13 +159,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
}
- /* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_stars_feedback) {
- if (cell_is_active_stars(ci, e)) {
- scheduler_activate(s, t);
- cell_activate_subcell_stars_tasks(ci, NULL, s);
- }
+ if (cell_is_active_stars(ci, e)) scheduler_activate(s, t);
}
/* Activate the gravity drift */
@@ -258,14 +251,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
/* Stars density */
- if (t_subtype == task_subtype_stars_density &&
- ((ci_active_stars && ci->nodeID == engine_rank) ||
- (cj_active_stars && cj->nodeID == engine_rank))) {
-
- // MATTHIEU: The logic here can be improved.
- // If ci is active for stars but not cj, then we can only drift the
- // stars in ci and parts in cj. (and vice-versa). The same logic can be
- // applied in cell_unskip_stars().
+ else if ((t_subtype == task_subtype_stars_density) &&
+ (ci_active_stars || cj_active_stars) &&
+ (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
scheduler_activate(s, t);
@@ -273,59 +261,65 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (t_type == task_type_pair) {
/* Do ci */
- /* Store some values. */
- atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
- atomic_or(&ci->stars.requires_sorts, 1 << t->flags);
-
- cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
- ci->stars.dx_max_sort_old = ci->stars.dx_max_sort;
+ if (ci_active_stars) {
- /* Activate the hydro drift tasks. */
- if (ci_nodeID == nodeID) cell_activate_drift_spart(ci, s);
+ /* stars for ci */
+ atomic_or(&ci->stars.requires_sorts, 1 << t->flags);
+ ci->stars.dx_max_sort_old = ci->stars.dx_max_sort;
- if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
+ /* hydro for cj */
+ atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
+ cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
- /* Check the sorts and activate them if needed. */
- cell_activate_hydro_sorts(cj, t->flags, s);
+ /* Activate the drift tasks. */
+ if (ci_nodeID == nodeID) cell_activate_drift_spart(ci, s);
+ if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
- cell_activate_stars_sorts(ci, t->flags, s);
+ /* Check the sorts and activate them if needed. */
+ cell_activate_hydro_sorts(cj, t->flags, s);
+ cell_activate_stars_sorts(ci, t->flags, s);
+ }
/* Do cj */
- /* Store some values. */
- atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
- atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
+ if (cj_active_stars) {
- ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
- cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
+ /* hydro for ci */
+ atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
+ ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
- /* Activate the hydro drift tasks. */
- if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
+ /* stars for cj */
+ atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
+ cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
- if (cj_nodeID == nodeID) cell_activate_drift_spart(cj, s);
+ /* Activate the drift tasks. */
+ if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
+ if (cj_nodeID == nodeID) cell_activate_drift_spart(cj, 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);
+ /* 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) {
- cell_activate_subcell_stars_tasks(t->ci, t->cj, s);
+ else if (t_type == task_type_sub_pair &&
+ t_subtype == task_subtype_stars_density) {
+ cell_activate_subcell_stars_tasks(ci, cj, s);
}
}
/* Stars feedback */
- if (t_subtype == task_subtype_stars_feedback &&
- ((ci_active_stars && ci->nodeID == engine_rank) ||
- (cj_active_stars && cj->nodeID == engine_rank))) {
+ else if ((t_subtype == task_subtype_stars_feedback) &&
+ ((ci_active_stars && ci_nodeID == nodeID) ||
+ (cj_active_stars && cj_nodeID == nodeID))) {
scheduler_activate(s, t);
}
/* Gravity */
- if ((t_subtype == task_subtype_grav) &&
- ((ci_active_gravity && ci_nodeID == nodeID) ||
- (cj_active_gravity && cj_nodeID == nodeID))) {
+ else if ((t_subtype == task_subtype_grav) &&
+ ((ci_active_gravity && ci_nodeID == nodeID) ||
+ (cj_active_gravity && cj_nodeID == nodeID))) {
scheduler_activate(s, t);
@@ -368,7 +362,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Is the foreign cell active and will need stuff from us? */
if (ci_active_hydro) {
-
struct link *l =
scheduler_activate_send(s, cj->mpi.hydro.send_xv, ci_nodeID);
@@ -396,6 +389,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* If the local cell is active, receive data from the foreign cell. */
if (ci_active_hydro) {
+
scheduler_activate(s, cj->mpi.hydro.recv_xv);
if (cj_active_hydro) {
scheduler_activate(s, cj->mpi.hydro.recv_rho);
@@ -439,16 +433,78 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
/* Only interested in stars_density tasks as of here. */
- if (t->subtype == task_subtype_stars_density) {
+ else if (t->subtype == task_subtype_stars_density) {
/* Too much particle movement? */
if (cell_need_rebuild_for_stars_pair(ci, cj)) *rebuild_space = 1;
+ if (cell_need_rebuild_for_stars_pair(cj, ci)) *rebuild_space = 1;
+
+#ifdef WITH_MPI
+ /* Activate the send/recv tasks. */
+ if (ci_nodeID != nodeID) {
+
+ if (cj_active_stars) {
+ scheduler_activate(s, ci->mpi.hydro.recv_xv);
+ scheduler_activate(s, ci->mpi.hydro.recv_rho);
- // LOIC: Need implementing MPI case
+ /* If the local cell is active, more stuff will be needed. */
+ scheduler_activate_send(s, cj->mpi.stars.send, ci_nodeID);
+ cell_activate_drift_spart(cj, s);
+
+ /* If the local cell is active, send its ti_end values. */
+ scheduler_activate_send(s, cj->mpi.send_ti, ci_nodeID);
+ }
+
+ if (ci_active_stars) {
+ scheduler_activate(s, ci->mpi.stars.recv);
+
+ /* If the foreign cell is active, we want its ti_end values. */
+ scheduler_activate(s, ci->mpi.recv_ti);
+
+ /* Is the foreign cell active and will need stuff from us? */
+ scheduler_activate_send(s, cj->mpi.hydro.send_xv, ci_nodeID);
+ scheduler_activate_send(s, cj->mpi.hydro.send_rho, ci_nodeID);
+
+ /* Drift the cell which will be sent; note that not all sent
+ particles will be drifted, only those that are needed. */
+ cell_activate_drift_part(cj, s);
+ }
+
+ } else if (cj_nodeID != nodeID) {
+
+ /* If the local cell is active, receive data from the foreign cell. */
+ if (ci_active_stars) {
+ scheduler_activate(s, cj->mpi.hydro.recv_xv);
+ scheduler_activate(s, cj->mpi.hydro.recv_rho);
+
+ /* If the local cell is active, more stuff will be needed. */
+ scheduler_activate_send(s, ci->mpi.stars.send, cj_nodeID);
+ cell_activate_drift_spart(ci, s);
+
+ /* If the local cell is active, send its ti_end values. */
+ scheduler_activate_send(s, ci->mpi.send_ti, cj_nodeID);
+ }
+
+ if (cj_active_stars) {
+ scheduler_activate(s, cj->mpi.stars.recv);
+
+ /* If the foreign cell is active, we want its ti_end values. */
+ scheduler_activate(s, cj->mpi.recv_ti);
+
+ /* Is the foreign cell active and will need stuff from us? */
+ scheduler_activate_send(s, ci->mpi.hydro.send_xv, cj_nodeID);
+ scheduler_activate_send(s, ci->mpi.hydro.send_rho, cj_nodeID);
+
+ /* Drift the cell which will be sent; note that not all sent
+ particles will be drifted, only those that are needed. */
+ cell_activate_drift_part(ci, s);
+ }
+ }
+#endif
}
/* Only interested in gravity tasks as of here. */
- if (t_subtype == task_subtype_grav) {
+ else if (t_subtype == task_subtype_grav) {
#ifdef WITH_MPI
/* Activate the send/recv tasks. */
@@ -504,11 +560,16 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
}
- /* End force ? */
- else if (t_type == task_type_end_force) {
+ /* End force for hydro ? */
+ else if (t_type == task_type_end_hydro_force) {
- if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e))
- scheduler_activate(s, t);
+ if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
+ }
+
+ /* End force for gravity ? */
+ else if (t_type == task_type_end_grav_force) {
+
+ if (cell_is_active_gravity(t->ci, e)) scheduler_activate(s, t);
}
/* Kick ? */
@@ -563,9 +624,12 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
/* Star ghost tasks ? */
- else if (t_type == task_type_stars_ghost ||
- t_type == task_type_stars_ghost_in ||
- t_type == task_type_stars_ghost_out) {
+ else if (t_type == task_type_stars_ghost) {
+ if (cell_is_active_stars(t->ci, e)) scheduler_activate(s, t);
+ }
+
+ /* Feedback implicit tasks? */
+ else if (t_type == task_type_stars_in || t_type == task_type_stars_out) {
if (cell_is_active_stars(t->ci, e)) scheduler_activate(s, t);
}
diff --git a/src/partition.c b/src/partition.c
index 1dde46dcd4ad2eb397bd26905fd8b57240088857..606f64e4c2f1c057520ed7ff893db10905f5aa8e 100644
--- a/src/partition.c
+++ b/src/partition.c
@@ -1244,9 +1244,10 @@ static void partition_gather_weights(void *map_data, int num_elements,
if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
t->type == task_type_ghost || t->type == task_type_extra_ghost ||
t->type == task_type_kick1 || t->type == task_type_kick2 ||
- t->type == task_type_end_force || t->type == task_type_cooling ||
- t->type == task_type_timestep || t->type == task_type_init_grav ||
- t->type == task_type_grav_down ||
+ t->type == task_type_end_hydro_force ||
+ t->type == task_type_end_grav_force || t->type == task_type_cooling ||
+ t->type == task_type_star_formation || t->type == task_type_timestep ||
+ t->type == task_type_init_grav || t->type == task_type_grav_down ||
t->type == task_type_grav_long_range) {
/* Particle updates add only to vertex weight. */
@@ -2149,9 +2150,10 @@ static void check_weights(struct task *tasks, int nr_tasks,
if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
t->type == task_type_ghost || t->type == task_type_extra_ghost ||
t->type == task_type_kick1 || t->type == task_type_kick2 ||
- t->type == task_type_end_force || t->type == task_type_cooling ||
- t->type == task_type_timestep || t->type == task_type_init_grav ||
- t->type == task_type_grav_down ||
+ t->type == task_type_end_hydro_force ||
+ t->type == task_type_end_grav_force || t->type == task_type_cooling ||
+ t->type == task_type_star_formation || t->type == task_type_timestep ||
+ t->type == task_type_init_grav || t->type == task_type_grav_down ||
t->type == task_type_grav_long_range) {
/* Particle updates add only to vertex weight. */
diff --git a/src/runner.c b/src/runner.c
index 5548bef3dd280f4f407f08eaff909918c4d36bcd..cfe3398684bef7b5a04cbfcefce1ebfe7729e24b 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -109,7 +109,9 @@
/* Import the stars density loop functions. */
#define FUNCTION density
+#define UPDATE_STARS 1
#include "runner_doiact_stars.h"
+#undef UPDATE_STARS
#undef FUNCTION
/* Import the stars feedback loop functions. */
@@ -130,11 +132,12 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
struct spart *restrict sparts = c->stars.parts;
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
- const struct stars_props *stars_properties = e->stars_properties;
- const float stars_h_max = stars_properties->h_max;
- const float eps = stars_properties->h_tolerance;
- const float stars_eta_dim = pow_dimension(stars_properties->eta_neighbours);
- const int max_smoothing_iter = stars_properties->max_smoothing_iterations;
+ const float stars_h_max = e->hydro_properties->h_max;
+ const float stars_h_min = e->hydro_properties->h_min;
+ const float eps = e->stars_properties->h_tolerance;
+ const float stars_eta_dim =
+ pow_dimension(e->stars_properties->eta_neighbours);
+ const int max_smoothing_iter = e->stars_properties->max_smoothing_iterations;
int redo = 0, scount = 0;
TIMER_TIC;
@@ -150,11 +153,23 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* Init the list of active particles that have to be updated. */
int *sid = NULL;
+ float *h_0 = NULL;
+ float *left = NULL;
+ float *right = NULL;
if ((sid = (int *)malloc(sizeof(int) * c->stars.count)) == NULL)
error("Can't allocate memory for sid.");
+ if ((h_0 = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
+ error("Can't allocate memory for h_0.");
+ if ((left = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
+ error("Can't allocate memory for left.");
+ if ((right = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
+ error("Can't allocate memory for right.");
for (int k = 0; k < c->stars.count; k++)
if (spart_is_active(&sparts[k], e)) {
sid[scount] = k;
+ h_0[scount] = sparts[k].h;
+ left[scount] = 0.f;
+ right[scount] = stars_h_max;
++scount;
}
@@ -178,9 +193,11 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
#endif
/* Get some useful values */
+ const float h_init = h_0[i];
const float h_old = sp->h;
const float h_old_dim = pow_dimension(h_old);
const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);
+
float h_new;
int has_no_neighbours = 0;
@@ -191,6 +208,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* Double h and try again */
h_new = 2.f * h_old;
+
} else {
/* Finish the density calculation */
@@ -204,8 +222,45 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
sp->density.wcount_dh * h_old_dim +
hydro_dimension * sp->density.wcount * h_old_dim_minus_one;
+ /* Improve the bisection bounds */
+ if (n_sum < n_target)
+ left[i] = max(left[i], h_old);
+ else if (n_sum > n_target)
+ right[i] = min(right[i], h_old);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check the validity of the left and right bounds */
+ if (left[i] > right[i])
+ error("Invalid left (%e) and right (%e)", left[i], right[i]);
+#endif
+
+ /* Skip if h is already h_max and we don't have enough neighbours */
+ /* Same if we are below h_min */
+ if (((sp->h >= stars_h_max) && (f < 0.f)) ||
+ ((sp->h <= stars_h_min) && (f > 0.f))) {
+
+ stars_reset_feedback(sp);
+
+ /* Ok, we are done with this particle */
+ continue;
+ }
+
+ /* Normal case: Use Newton-Raphson to get a better value of h */
+
/* Avoid floating point exception from f_prime = 0 */
h_new = h_old - f / (f_prime + FLT_MIN);
+
+ /* Be verbose about the particles that struggle to converge */
+ if (num_reruns > max_smoothing_iter - 10) {
+
+ message(
+ "Smoothing length convergence problem: iter=%d p->id=%lld "
+ "h_init=%12.8e h_old=%12.8e h_new=%12.8e f=%f f_prime=%f "
+ "n_sum=%12.8e n_target=%12.8e left=%12.8e right=%12.8e",
+ num_reruns, sp->id, h_init, h_old, h_new, f, f_prime, n_sum,
+ n_target, left[i], right[i]);
+ }
+
#ifdef SWIFT_DEBUG_CHECKS
if ((f > 0.f && h_new > h_old) || (f < 0.f && h_new < h_old))
error(
@@ -215,19 +270,39 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
h_new = min(h_new, 2.f * h_old);
h_new = max(h_new, 0.5f * h_old);
+
+ /* Verify that we are actually progrssing towards the answer */
+ h_new = max(h_new, left[i]);
+ h_new = min(h_new, right[i]);
}
/* Check whether the particle has an inappropriate smoothing length */
if (fabsf(h_new - h_old) > eps * h_old) {
/* Ok, correct then */
- sp->h = h_new;
+
+ /* Case where we have been oscillating around the solution */
+ if ((h_new == left[i] && h_old == right[i]) ||
+ (h_old == left[i] && h_new == right[i])) {
+
+ /* Bissect the remaining interval */
+ sp->h = pow_inv_dimension(
+ 0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));
+
+ } else {
+
+ /* Normal case */
+ sp->h = h_new;
+ }
/* If below the absolute maximum, try again */
- if (sp->h < stars_h_max) {
+ if (sp->h < stars_h_max && sp->h > stars_h_min) {
/* Flag for another round of fun */
sid[redo] = sid[i];
+ h_0[redo] = h_0[i];
+ left[redo] = left[i];
+ right[redo] = right[i];
redo += 1;
/* Re-initialise everything */
@@ -236,7 +311,12 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* Off we go ! */
continue;
- } else {
+ } else if (sp->h <= stars_h_min) {
+
+ /* Ok, this particle is a lost cause... */
+ sp->h = stars_h_min;
+
+ } else if (sp->h >= stars_h_max) {
/* Ok, this particle is a lost cause... */
sp->h = stars_h_max;
@@ -245,13 +325,19 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
if (has_no_neighbours) {
stars_spart_has_no_neighbours(sp, cosmo);
}
+
+ } else {
+ error(
+ "Fundamental problem with the smoothing length iteration "
+ "logic.");
}
}
/* We now have a particle whose smoothing length has converged */
+ stars_reset_feedback(sp);
/* Compute the stellar evolution */
- stars_evolve_spart(sp, stars_properties, cosmo);
+ stars_evolve_spart(sp, e->stars_properties, cosmo);
}
/* We now need to treat the particles whose smoothing length had not
@@ -315,7 +401,10 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
}
/* Be clean */
+ free(left);
+ free(right);
free(sid);
+ free(h_0);
}
if (timer) TIMER_TOC(timer_dostars_ghost);
@@ -488,6 +577,7 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
struct cooling_function_data *restrict cooling = e->cooling_func;
const double time_base = e->time_base;
const integertime_t ti_current = e->ti_current;
+ const int current_stars_count = c->stars.count;
TIMER_TIC;
@@ -555,6 +645,13 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
} /* Loop over particles */
}
+ /* If we formed any stars, the star sorts are now invalid. We need to
+ * re-compute them. */
+ if ((c == c->hydro.super) && (current_stars_count != c->stars.count)) {
+ cell_clear_stars_sort_flags(c, /*is_super=*/1);
+ runner_do_stars_sort(r, c, 0x1FFF, /*cleanup=*/0, /*timer=*/0);
+ }
+
if (timer) TIMER_TOC(timer_do_star_formation);
}
@@ -917,8 +1014,9 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
if (flags == 0 && !c->stars.do_sub_sort) return;
/* Check that the particles have been moved to the current time */
- if (flags && !cell_are_spart_drifted(c, r->e))
+ if (flags && !cell_are_spart_drifted(c, r->e)) {
error("Sorting un-drifted cell c->nodeID=%d", c->nodeID);
+ }
#ifdef SWIFT_DEBUG_CHECKS
/* Make sure the sort flags are consistent (downward). */
@@ -954,10 +1052,10 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL && c->progeny[k]->stars.count > 0) {
/* Only propagate cleanup if the progeny is stale. */
- runner_do_stars_sort(r, c->progeny[k], flags,
- cleanup && (c->progeny[k]->stars.dx_max_sort_old >
- space_maxreldx * c->progeny[k]->dmin),
- 0);
+ const int cleanup_prog =
+ cleanup && (c->progeny[k]->stars.dx_max_sort_old >
+ space_maxreldx * c->progeny[k]->dmin);
+ runner_do_stars_sort(r, c->progeny[k], flags, cleanup_prog, 0);
dx_max_sort = max(dx_max_sort, c->progeny[k]->stars.dx_max_sort);
dx_max_sort_old =
max(dx_max_sort_old, c->progeny[k]->stars.dx_max_sort_old);
@@ -1324,8 +1422,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
hydro_dimension * p->density.wcount * h_old_dim_minus_one;
/* Improve the bisection bounds */
- if (n_sum < n_target) left[i] = max(left[i], h_old);
- if (n_sum > n_target) right[i] = min(right[i], h_old);
+ if (n_sum < n_target)
+ left[i] = max(left[i], h_old);
+ else if (n_sum > n_target)
+ right[i] = min(right[i], h_old);
#ifdef SWIFT_DEBUG_CHECKS
/* Check the validity of the left and right bounds */
@@ -1334,6 +1434,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
#endif
/* Skip if h is already h_max and we don't have enough neighbours */
+ /* Same if we are below h_min */
if (((p->h >= hydro_h_max) && (f < 0.f)) ||
((p->h <= hydro_h_min) && (f > 0.f))) {
@@ -1474,6 +1575,11 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
hydro_part_has_no_neighbours(p, xp, cosmo);
chemistry_part_has_no_neighbours(p, xp, chemistry, cosmo);
}
+
+ } else {
+ error(
+ "Fundamental problem with the smoothing length iteration "
+ "logic.");
}
}
@@ -1711,6 +1817,7 @@ void runner_do_unskip_mapper(void *map_data, int num_elements,
}
}
}
+
/**
* @brief Drift all part in a cell.
*
@@ -1743,6 +1850,21 @@ void runner_do_drift_gpart(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_drift_gpart);
}
+/**
+ * @brief Drift all spart in a cell.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_drift_spart(struct runner *r, struct cell *c, int timer) {
+
+ TIMER_TIC;
+
+ cell_drift_spart(c, r->e, 0);
+
+ if (timer) TIMER_TOC(timer_drift_spart);
+}
/**
* @brief Perform the first half-kick on all the active particles in a cell.
*
@@ -1770,7 +1892,9 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
/* Anything to do here? */
- if (!cell_is_starting_hydro(c, e) && !cell_is_starting_gravity(c, e)) return;
+ if (!cell_is_starting_hydro(c, e) && !cell_is_starting_gravity(c, e) &&
+ !cell_is_starting_stars(c, e))
+ return;
/* Recurse? */
if (c->split) {
@@ -1952,7 +2076,9 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
/* Anything to do here? */
- if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e)) return;
+ if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e) &&
+ !cell_is_active_stars(c, e))
+ return;
/* Recurse? */
if (c->split) {
@@ -2146,7 +2272,8 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
/* Anything to do here? */
- if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e)) {
+ if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e) &&
+ !cell_is_active_stars(c, e)) {
c->hydro.updated = 0;
c->grav.updated = 0;
c->stars.updated = 0;
@@ -2159,7 +2286,8 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
ti_hydro_beg_max = 0;
integertime_t ti_gravity_end_min = max_nr_timesteps, ti_gravity_end_max = 0,
ti_gravity_beg_max = 0;
- integertime_t ti_stars_end_min = max_nr_timesteps;
+ integertime_t ti_stars_end_min = max_nr_timesteps, ti_stars_end_max = 0,
+ ti_stars_beg_max = 0;
/* No children? */
if (!c->split) {
@@ -2337,13 +2465,13 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
s_updated++;
g_updated++;
- /* What is the next sync-point ? */
+ ti_stars_end_min = min(ti_current + ti_new_step, ti_stars_end_min);
+ ti_stars_end_max = max(ti_current + ti_new_step, ti_stars_end_max);
ti_gravity_end_min = min(ti_current + ti_new_step, ti_gravity_end_min);
ti_gravity_end_max = max(ti_current + ti_new_step, ti_gravity_end_max);
- ti_stars_end_min = min(ti_current + ti_new_step, ti_stars_end_min);
-
/* What is the next starting point for this cell ? */
+ ti_stars_beg_max = max(ti_current, ti_stars_beg_max);
ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
/* star particle is inactive but not inhibited */
@@ -2355,23 +2483,23 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
const integertime_t ti_end =
get_integer_time_end(ti_current, sp->time_bin);
- /* What is the next sync-point ? */
- ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
- ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
-
- ti_stars_end_min = min(ti_end, ti_stars_end_min);
-
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, sp->time_bin);
+ ti_stars_end_min = min(ti_end, ti_stars_end_min);
+ ti_stars_end_max = max(ti_end, ti_stars_end_max);
+ ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
+ ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
+
/* What is the next starting point for this cell ? */
+ ti_stars_beg_max = max(ti_beg, ti_stars_beg_max);
ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
}
}
} else {
/* Loop over the progeny. */
- for (int k = 0; k < 8; k++)
+ for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
struct cell *restrict cp = c->progeny[k];
@@ -2392,7 +2520,10 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
ti_gravity_end_max = max(cp->grav.ti_end_max, ti_gravity_end_max);
ti_gravity_beg_max = max(cp->grav.ti_beg_max, ti_gravity_beg_max);
ti_stars_end_min = min(cp->stars.ti_end_min, ti_stars_end_min);
+ ti_stars_end_max = max(cp->grav.ti_end_max, ti_stars_end_max);
+ ti_stars_beg_max = max(cp->grav.ti_beg_max, ti_stars_beg_max);
}
+ }
}
/* Store the values. */
@@ -2409,6 +2540,8 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
c->grav.ti_end_max = ti_gravity_end_max;
c->grav.ti_beg_max = ti_gravity_beg_max;
c->stars.ti_end_min = ti_stars_end_min;
+ c->stars.ti_end_max = ti_stars_end_max;
+ c->stars.ti_beg_max = ti_stars_beg_max;
#ifdef SWIFT_DEBUG_CHECKS
if (c->hydro.ti_end_min == e->ti_current &&
@@ -2564,46 +2697,32 @@ void runner_do_limiter(struct runner *r, struct cell *c, int force, int timer) {
}
/**
- * @brief End the force calculation of all active particles in a cell
+ * @brief End the hydro force calculation of all active particles in a cell
* by multiplying the acccelerations by the relevant constants
*
* @param r The #runner thread.
* @param c The #cell.
* @param timer Are we timing this ?
*/
-void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
+void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
- const struct space *s = e->s;
- const struct cosmology *cosmo = e->cosmology;
- const int count = c->hydro.count;
- const int gcount = c->grav.count;
- const int scount = c->stars.count;
- struct part *restrict parts = c->hydro.parts;
- struct gpart *restrict gparts = c->grav.parts;
- struct spart *restrict sparts = c->stars.parts;
- const int periodic = s->periodic;
- const float G_newton = e->physical_constants->const_newton_G;
TIMER_TIC;
- /* Potential normalisation in the case of periodic gravity */
- float potential_normalisation = 0.;
- if (periodic && (e->policy & engine_policy_self_gravity)) {
- const double volume = s->dim[0] * s->dim[1] * s->dim[2];
- const double r_s = e->mesh->r_s;
- potential_normalisation = 4. * M_PI * e->total_mass * r_s * r_s / volume;
- }
-
/* Anything to do here? */
- if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e)) return;
+ if (!cell_is_active_hydro(c, e)) return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) runner_do_end_force(r, c->progeny[k], 0);
+ if (c->progeny[k] != NULL) runner_do_end_hydro_force(r, c->progeny[k], 0);
} else {
+ const struct cosmology *cosmo = e->cosmology;
+ const int count = c->hydro.count;
+ struct part *restrict parts = c->hydro.parts;
+
/* Loop over the gas particles in this cell. */
for (int k = 0; k < count; k++) {
@@ -2616,6 +2735,48 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
hydro_end_force(p, cosmo);
}
}
+ }
+
+ if (timer) TIMER_TOC(timer_end_hydro_force);
+}
+
+/**
+ * @brief End the gravity force calculation of all active particles in a cell
+ * by multiplying the acccelerations by the relevant constants
+ *
+ * @param r The #runner thread.
+ * @param c The #cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_end_grav_force(struct runner *r, struct cell *c, int timer) {
+
+ const struct engine *e = r->e;
+
+ TIMER_TIC;
+
+ /* Anything to do here? */
+ if (!cell_is_active_gravity(c, e)) return;
+
+ /* Recurse? */
+ if (c->split) {
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) runner_do_end_grav_force(r, c->progeny[k], 0);
+ } else {
+
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const float G_newton = e->physical_constants->const_newton_G;
+
+ /* Potential normalisation in the case of periodic gravity */
+ float potential_normalisation = 0.;
+ if (periodic && (e->policy & engine_policy_self_gravity)) {
+ const double volume = s->dim[0] * s->dim[1] * s->dim[2];
+ const double r_s = e->mesh->r_s;
+ potential_normalisation = 4. * M_PI * e->total_mass * r_s * r_s / volume;
+ }
+
+ const int gcount = c->grav.count;
+ struct gpart *restrict gparts = c->grav.parts;
/* Loop over the g-particles in this cell. */
for (int k = 0; k < gcount; k++) {
@@ -2693,21 +2854,8 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
#endif
}
}
-
- /* Loop over the stars particles in this cell. */
- for (int k = 0; k < scount; k++) {
-
- /* Get a handle on the spart. */
- struct spart *restrict sp = &sparts[k];
- if (spart_is_active(sp, e)) {
-
- /* Finish the force loop */
- stars_end_force(sp);
- }
- }
}
-
- if (timer) TIMER_TOC(timer_endforce);
+ if (timer) TIMER_TOC(timer_end_grav_force);
}
/**
@@ -2720,7 +2868,6 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
*/
void runner_do_recv_part(struct runner *r, struct cell *c, int clear_sorts,
int timer) {
-
#ifdef WITH_MPI
const struct part *restrict parts = c->hydro.parts;
@@ -2872,69 +3019,79 @@ void runner_do_recv_gpart(struct runner *r, struct cell *c, int timer) {
*
* @param r The runner thread.
* @param c The cell.
+ * @param clear_sorts Should we clear the sort flag and hence trigger a sort ?
* @param timer Are we timing this ?
*/
-void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
+void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
+ int timer) {
#ifdef WITH_MPI
- const struct spart *restrict sparts = c->stars.parts;
+ struct spart *restrict sparts = c->stars.parts;
const size_t nr_sparts = c->stars.count;
const integertime_t ti_current = r->e->ti_current;
- error("Need to add h_max computation");
-
TIMER_TIC;
- integertime_t ti_gravity_end_min = max_nr_timesteps;
- integertime_t ti_gravity_end_max = 0;
+ integertime_t ti_stars_end_min = max_nr_timesteps;
+ integertime_t ti_stars_end_max = 0;
timebin_t time_bin_min = num_time_bins;
timebin_t time_bin_max = 0;
+ float h_max = 0.f;
#ifdef SWIFT_DEBUG_CHECKS
if (c->nodeID == engine_rank) error("Updating a local cell!");
#endif
+ /* Clear this cell's sorted mask. */
+ if (clear_sorts) c->stars.sorted = 0;
+
/* If this cell is a leaf, collect the particle data. */
if (!c->split) {
/* Collect everything... */
for (size_t k = 0; k < nr_sparts; k++) {
+#ifdef DEBUG_INTERACTIONS_STARS
+ sparts[k].num_ngb_force = 0;
+#endif
if (sparts[k].time_bin == time_bin_inhibited) continue;
time_bin_min = min(time_bin_min, sparts[k].time_bin);
time_bin_max = max(time_bin_max, sparts[k].time_bin);
+ h_max = max(h_max, sparts[k].h);
}
/* Convert into a time */
- ti_gravity_end_min = get_integer_time_end(ti_current, time_bin_min);
- ti_gravity_end_max = get_integer_time_end(ti_current, time_bin_max);
+ ti_stars_end_min = get_integer_time_end(ti_current, time_bin_min);
+ ti_stars_end_max = get_integer_time_end(ti_current, time_bin_max);
}
/* Otherwise, recurse and collect. */
else {
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL && c->progeny[k]->stars.count > 0) {
- runner_do_recv_spart(r, c->progeny[k], 0);
- ti_gravity_end_min =
- min(ti_gravity_end_min, c->progeny[k]->grav.ti_end_min);
- ti_gravity_end_max =
- max(ti_gravity_end_max, c->progeny[k]->grav.ti_end_max);
+ runner_do_recv_spart(r, c->progeny[k], clear_sorts, 0);
+ ti_stars_end_min =
+ min(ti_stars_end_min, c->progeny[k]->stars.ti_end_min);
+ ti_stars_end_max =
+ max(ti_stars_end_max, c->progeny[k]->grav.ti_end_max);
+ h_max = max(h_max, c->progeny[k]->stars.h_max);
}
}
}
#ifdef SWIFT_DEBUG_CHECKS
- if (ti_gravity_end_min < ti_current)
+ if (ti_stars_end_min < ti_current)
error(
"Received a cell at an incorrect time c->ti_end_min=%lld, "
"e->ti_current=%lld.",
- ti_gravity_end_min, ti_current);
+ ti_stars_end_min, ti_current);
#endif
/* ... and store. */
// c->grav.ti_end_min = ti_gravity_end_min;
// c->grav.ti_end_max = ti_gravity_end_max;
- c->grav.ti_old_part = ti_current;
+ c->stars.ti_old_part = ti_current;
+ c->stars.h_max = h_max;
if (timer) TIMER_TOC(timer_dorecv_spart);
@@ -3022,9 +3179,9 @@ void *runner_main(void *data) {
else if (t->subtype == task_subtype_external_grav)
runner_do_grav_external(r, ci, 1);
else if (t->subtype == task_subtype_stars_density)
- runner_doself_stars_density(r, ci, 1);
+ runner_doself_branch_stars_density(r, ci);
else if (t->subtype == task_subtype_stars_feedback)
- runner_doself_stars_feedback(r, ci, 1);
+ runner_doself_branch_stars_feedback(r, ci);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
@@ -3043,9 +3200,9 @@ void *runner_main(void *data) {
else if (t->subtype == task_subtype_grav)
runner_dopair_recursive_grav(r, ci, cj, 1);
else if (t->subtype == task_subtype_stars_density)
- runner_dopair_stars_density(r, ci, cj, 1);
+ runner_dopair_branch_stars_density(r, ci, cj);
else if (t->subtype == task_subtype_stars_feedback)
- runner_dopair_stars_feedback(r, ci, cj, 1);
+ runner_dopair_branch_stars_feedback(r, ci, cj);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
@@ -3121,6 +3278,9 @@ void *runner_main(void *data) {
case task_type_drift_part:
runner_do_drift_part(r, ci, 1);
break;
+ case task_type_drift_spart:
+ runner_do_drift_spart(r, ci, 1);
+ break;
case task_type_drift_gpart:
runner_do_drift_gpart(r, ci, 1);
break;
@@ -3130,8 +3290,11 @@ void *runner_main(void *data) {
case task_type_kick2:
runner_do_kick2(r, ci, 1);
break;
- case task_type_end_force:
- runner_do_end_force(r, ci, 1);
+ case task_type_end_hydro_force:
+ runner_do_end_hydro_force(r, ci, 1);
+ break;
+ case task_type_end_grav_force:
+ runner_do_end_grav_force(r, ci, 1);
break;
case task_type_logger:
runner_do_logger(r, ci, 1);
@@ -3163,7 +3326,7 @@ void *runner_main(void *data) {
} else if (t->subtype == task_subtype_gpart) {
runner_do_recv_gpart(r, ci, 1);
} else if (t->subtype == task_subtype_spart) {
- runner_do_recv_spart(r, ci, 1);
+ runner_do_recv_spart(r, ci, 1, 1);
} else if (t->subtype == task_subtype_multipole) {
cell_unpack_multipoles(ci, (struct gravity_tensors *)t->buff);
free(t->buff);
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index 861798b70b8ba90b9267375253bd8570baec3e9a..383fa64c56f6e6e7a134452f9eef8f715171be0c 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -803,6 +803,9 @@ void DOPAIR_SUBSET_BRANCH(struct runner *r, struct cell *restrict ci,
const struct engine *e = r->e;
+ /* Anything to do here? */
+ if (cj->hydro.count == 0) return;
+
/* Get the relative distance between the pairs, wrapping. */
double shift[3] = {0.0, 0.0, 0.0};
for (int k = 0; k < 3; k++) {
@@ -1187,6 +1190,9 @@ void DOPAIR1_BRANCH(struct runner *r, struct cell *ci, struct cell *cj) {
const struct engine *restrict e = r->e;
+ /* Anything to do here? */
+ if (ci->hydro.count == 0 || cj->hydro.count == 0) return;
+
/* Anything to do here? */
if (!cell_is_active_hydro(ci, e) && !cell_is_active_hydro(cj, e)) return;
@@ -1725,6 +1731,9 @@ void DOPAIR2_BRANCH(struct runner *r, struct cell *ci, struct cell *cj) {
const struct engine *restrict e = r->e;
+ /* Anything to do here? */
+ if (ci->hydro.count == 0 || cj->hydro.count == 0) return;
+
/* Anything to do here? */
if (!cell_is_active_hydro(ci, e) && !cell_is_active_hydro(cj, e)) return;
@@ -1971,6 +1980,9 @@ void DOSELF1_BRANCH(struct runner *r, struct cell *c) {
const struct engine *restrict e = r->e;
+ /* Anything to do here? */
+ if (c->hydro.count == 0) return;
+
/* Anything to do here? */
if (!cell_is_active_hydro(c, e)) return;
@@ -2145,6 +2157,9 @@ void DOSELF2_BRANCH(struct runner *r, struct cell *c) {
const struct engine *restrict e = r->e;
+ /* Anything to do here? */
+ if (c->hydro.count == 0) return;
+
/* Anything to do here? */
if (!cell_is_active_hydro(c, e)) return;
diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h
index e816d80399a0fef85645c914168dd4038f55988c..f208f14ac98a31a55df6741a2cc5f9cb0a829762 100644
--- a/src/runner_doiact_stars.h
+++ b/src/runner_doiact_stars.h
@@ -74,6 +74,11 @@
* @param timer 1 if the time is to be recorded.
*/
void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->nodeID != engine_rank) error("Should be run on a different node");
+#endif
+
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
@@ -95,6 +100,8 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
/* Get a hold of the ith spart in ci. */
struct spart *restrict si = &sparts[sid];
+ if (!spart_is_active(si, e) || spart_is_inhibited(si, e)) continue;
+
const float hi = si->h;
const float hig2 = hi * hi * kernel_gamma2;
const float six[3] = {(float)(si->x[0] - c->loc[0]),
@@ -138,6 +145,14 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
struct cell *restrict cj) {
+#ifdef SWIFT_DEBUG_CHECKS
+#ifdef UPDATE_STARS
+ if (ci->nodeID != engine_rank) error("Should be run on a different node");
+#else
+ if (cj->nodeID != engine_rank) error("Should be run on a different node");
+#endif
+#endif
+
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
@@ -167,6 +182,7 @@ void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
/* Get a hold of the ith spart in ci. */
struct spart *restrict si = &sparts_i[sid];
+ if (!spart_is_active(si, e) || spart_is_inhibited(si, e)) continue;
const float hi = si->h;
const float hig2 = hi * hi * kernel_gamma2;
const float six[3] = {(float)(si->x[0] - (cj->loc[0] + shift[0])),
@@ -202,9 +218,17 @@ void DO_NONSYM_PAIR1_STARS(struct runner *r, struct cell *restrict ci,
void DOPAIR1_STARS(struct runner *r, struct cell *restrict ci,
struct cell *restrict cj, int timer) {
- if (ci->stars.count != 0 && cj->hydro.count != 0)
+#ifdef UPDATE_STARS
+ const int ci_local = ci->nodeID == engine_rank;
+ const int cj_local = cj->nodeID == engine_rank;
+#else
+ /* here we are updating the hydro -> switch ci, cj */
+ const int ci_local = cj->nodeID == engine_rank;
+ const int cj_local = ci->nodeID == engine_rank;
+#endif
+ if (ci_local && ci->stars.count != 0 && cj->hydro.count != 0)
DO_NONSYM_PAIR1_STARS(r, ci, cj);
- if (cj->stars.count != 0 && ci->hydro.count != 0)
+ if (cj_local && cj->stars.count != 0 && ci->hydro.count != 0)
DO_NONSYM_PAIR1_STARS(r, cj, ci);
}
@@ -227,6 +251,9 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
int scount, struct cell *restrict cj,
const double *shift) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci->nodeID != engine_rank) error("Should be run on a different node");
+#endif
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
@@ -293,6 +320,10 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
struct spart *restrict sparts, int *restrict ind,
int scount) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci->nodeID != engine_rank) error("Should be run on a different node");
+#endif
+
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
@@ -446,8 +477,8 @@ void DOSUB_SUBSET_STARS(struct runner *r, struct cell *ci, struct spart *sparts,
else {
/* Recurse? */
- if (cell_can_recurse_in_pair_stars_task(ci) &&
- cell_can_recurse_in_pair_stars_task(cj)) {
+ if (cell_can_recurse_in_pair_stars_task(ci, cj) &&
+ cell_can_recurse_in_pair_stars_task(cj, ci)) {
/* Get the type of pair if not specified explicitly. */
double shift[3] = {0.0, 0.0, 0.0};
@@ -956,10 +987,13 @@ void DOSUB_SUBSET_STARS(struct runner *r, struct cell *ci, struct spart *sparts,
}
/* Otherwise, compute the pair directly. */
- else if (cell_is_active_stars(ci, e) || cell_is_active_stars(cj, e)) {
+ else if (cell_is_active_stars(ci, e)) {
/* Do any of the cells need to be drifted first? */
- if (!cell_are_part_drifted(cj, e)) error("Cell should be drifted!");
+ if (cell_is_active_stars(ci, e)) {
+ if (!cell_are_spart_drifted(ci, e)) error("Cell should be drifted!");
+ if (!cell_are_part_drifted(cj, e)) error("Cell should be drifted!");
+ }
DOPAIR1_SUBSET_BRANCH_STARS(r, ci, sparts, ind, scount, cj);
}
@@ -979,6 +1013,9 @@ void DOSELF1_BRANCH_STARS(struct runner *r, struct cell *c) {
const struct engine *restrict e = r->e;
+ /* Anything to do here? */
+ if (c->stars.count == 0) return;
+
/* Anything to do here? */
if (!cell_is_active_stars(c, e)) return;
@@ -995,6 +1032,8 @@ void DOSELF1_BRANCH_STARS(struct runner *r, struct cell *c) {
\
for (int pjd = 0; pjd < cj->TYPE.count; pjd++) { \
const struct PART *p = &cj->TYPE.parts[sort_j[pjd].i]; \
+ if (PART##_is_inhibited(p, e)) continue; \
+ \
const float d = p->x[0] * runner_shift[sid][0] + \
p->x[1] * runner_shift[sid][1] + \
p->x[2] * runner_shift[sid][2]; \
@@ -1007,9 +1046,12 @@ void DOSELF1_BRANCH_STARS(struct runner *r, struct cell *c) {
"cj->nodeID=%d " \
"ci->nodeID=%d d=%e sort_j[pjd].d=%e cj->" #TYPE \
".dx_max_sort=%e " \
- "cj->" #TYPE ".dx_max_sort_old=%e", \
+ "cj->" #TYPE \
+ ".dx_max_sort_old=%e, cellID=%i super->cellID=%i" \
+ "cj->depth=%d cj->maxdepth=%d", \
cj->nodeID, ci->nodeID, d, sort_j[pjd].d, cj->TYPE.dx_max_sort, \
- cj->TYPE.dx_max_sort_old); \
+ cj->TYPE.dx_max_sort_old, cj->cellID, cj->hydro.super->cellID, \
+ cj->depth, cj->maxdepth); \
} \
})
@@ -1026,18 +1068,29 @@ void DOSELF1_BRANCH_STARS(struct runner *r, struct cell *c) {
void DOPAIR1_BRANCH_STARS(struct runner *r, struct cell *ci, struct cell *cj) {
const struct engine *restrict e = r->e;
+
+ /* Get the sort ID. */
+ double shift[3] = {0.0, 0.0, 0.0};
+ const int sid = space_getsid(e->s, &ci, &cj, shift);
+
const int ci_active = cell_is_active_stars(ci, e);
const int cj_active = cell_is_active_stars(cj, e);
- const int do_ci = (ci->stars.count != 0 && cj->hydro.count != 0 && ci_active);
- const int do_cj = (cj->stars.count != 0 && ci->hydro.count != 0 && cj_active);
+#ifdef UPDATE_STARS
+ const int ci_local = ci->nodeID == engine_rank;
+ const int cj_local = cj->nodeID == engine_rank;
+#else
+ /* here we are updating the hydro -> switch ci, cj */
+ const int ci_local = cj->nodeID == engine_rank;
+ const int cj_local = ci->nodeID == engine_rank;
+#endif
+ const int do_ci =
+ (ci->stars.count != 0 && cj->hydro.count != 0 && ci_active && ci_local);
+ const int do_cj =
+ (cj->stars.count != 0 && ci->hydro.count != 0 && cj_active && cj_local);
/* Anything to do here? */
if (!do_ci && !do_cj) return;
- /* Get the sort ID. */
- double shift[3] = {0.0, 0.0, 0.0};
- const int sid = space_getsid(e->s, &ci, &cj, shift);
-
/* Check that cells are drifted. */
if (do_ci &&
(!cell_are_spart_drifted(ci, e) || !cell_are_part_drifted(cj, e)))
@@ -1099,19 +1152,19 @@ void DOSUB_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
const struct engine *e = r->e;
/* Should we even bother? */
- int should_do = ci->stars.count != 0 && cj->hydro.count != 0 &&
- cell_is_active_stars(ci, e);
- should_do |= cj->stars.count != 0 && ci->hydro.count != 0 &&
- cell_is_active_stars(cj, e);
- if (!should_do) return;
+ const int should_do_ci = ci->stars.count != 0 && cj->hydro.count != 0 &&
+ cell_is_active_stars(ci, e);
+ const int should_do_cj = cj->stars.count != 0 && ci->hydro.count != 0 &&
+ cell_is_active_stars(cj, e);
+ if (!should_do_ci && !should_do_cj) return;
/* Get the type of pair if not specified explicitly. */
double shift[3];
sid = space_getsid(s, &ci, &cj, shift);
/* Recurse? */
- if (cell_can_recurse_in_pair_stars_task(ci) &&
- cell_can_recurse_in_pair_stars_task(cj)) {
+ if (cell_can_recurse_in_pair_stars_task(ci, cj) &&
+ cell_can_recurse_in_pair_stars_task(cj, ci)) {
/* Different types of flags. */
switch (sid) {
@@ -1314,10 +1367,18 @@ void DOSUB_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
/* Otherwise, compute the pair directly. */
else {
+#ifdef UPDATE_STARS
+ const int ci_local = ci->nodeID == engine_rank;
+ const int cj_local = cj->nodeID == engine_rank;
+#else
+ /* here we are updating the hydro -> switch ci, cj */
+ const int ci_local = cj->nodeID == engine_rank;
+ const int cj_local = ci->nodeID == engine_rank;
+#endif
const int do_ci = ci->stars.count != 0 && cj->hydro.count != 0 &&
- cell_is_active_stars(ci, e);
+ cell_is_active_stars(ci, e) && ci_local;
const int do_cj = cj->stars.count != 0 && ci->hydro.count != 0 &&
- cell_is_active_stars(cj, e);
+ cell_is_active_stars(cj, e) && cj_local;
if (do_ci) {
@@ -1330,8 +1391,9 @@ void DOSUB_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
/* Do any of the cells need to be sorted first? */
if (!(ci->stars.sorted & (1 << sid)) ||
- ci->stars.dx_max_sort_old > ci->dmin * space_maxreldx)
+ ci->stars.dx_max_sort_old > ci->dmin * space_maxreldx) {
error("Interacting unsorted cell (sparts).");
+ }
if (!(cj->hydro.sorted & (1 << sid)) ||
cj->hydro.dx_max_sort_old > cj->dmin * space_maxreldx)
@@ -1349,12 +1411,14 @@ void DOSUB_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
/* Do any of the cells need to be sorted first? */
if (!(ci->hydro.sorted & (1 << sid)) ||
- ci->hydro.dx_max_sort_old > ci->dmin * space_maxreldx)
+ ci->hydro.dx_max_sort_old > ci->dmin * space_maxreldx) {
error("Interacting unsorted cell (parts).");
+ }
if (!(cj->stars.sorted & (1 << sid)) ||
- cj->stars.dx_max_sort_old > cj->dmin * space_maxreldx)
+ cj->stars.dx_max_sort_old > cj->dmin * space_maxreldx) {
error("Interacting unsorted cell (sparts).");
+ }
}
if (do_ci || do_cj) DOPAIR1_BRANCH_STARS(r, ci, cj);
@@ -1370,6 +1434,11 @@ void DOSUB_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
*/
void DOSUB_SELF1_STARS(struct runner *r, struct cell *ci, int gettimer) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci->nodeID != engine_rank)
+ error("This function should not be called on foreign cells");
+#endif
+
/* Should we even bother? */
if (ci->hydro.count == 0 || ci->stars.count == 0 ||
!cell_is_active_stars(ci, r->e))
diff --git a/src/scheduler.c b/src/scheduler.c
index 1083be1996f9ba3986a601ecbe47b1e6a224239a..229c94b2088814352487032816f610a1bdf74724 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -541,6 +541,11 @@ void scheduler_write_dependencies(struct scheduler *s, int verbose) {
*/
static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
+ /* Are we considering both stars and hydro when splitting? */
+ /* Note this is not very clean as the scheduler should not really
+ access the engine... */
+ const int with_feedback = (s->space->e->policy & engine_policy_feedback);
+
/* Iterate on this task until we're done with it. */
int redo = 1;
while (redo) {
@@ -548,9 +553,20 @@ static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
/* Reset the redo flag. */
redo = 0;
+ /* Is this a non-empty self-task? */
+ const int is_self =
+ (t->type == task_type_self) && (t->ci != NULL) &&
+ ((t->ci->hydro.count > 0) || (with_feedback && t->ci->stars.count > 0));
+
+ /* Is this a non-empty pair-task? */
+ const int is_pair =
+ (t->type == task_type_pair) && (t->ci != NULL) && (t->cj != NULL) &&
+ ((t->ci->hydro.count > 0) ||
+ (with_feedback && t->ci->stars.count > 0)) &&
+ ((t->cj->hydro.count > 0) || (with_feedback && t->cj->stars.count > 0));
+
/* Empty task? */
- if ((t->ci == NULL) || (t->type == task_type_pair && t->cj == NULL) ||
- t->ci->hydro.count == 0 || (t->cj != NULL && t->cj->hydro.count == 0)) {
+ if (!is_self && !is_pair) {
t->type = task_type_none;
t->subtype = task_subtype_none;
t->cj = NULL;
@@ -589,24 +605,46 @@ static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
int first_child = 0;
while (ci->progeny[first_child] == NULL) first_child++;
t->ci = ci->progeny[first_child];
- for (int k = first_child + 1; k < 8; k++)
- if (ci->progeny[k] != NULL && ci->progeny[k]->hydro.count)
+ for (int k = first_child + 1; k < 8; k++) {
+
+ /* Do we have a non-empty progenitor? */
+ if (ci->progeny[k] != NULL &&
+ (ci->progeny[k]->hydro.count ||
+ (with_feedback && ci->progeny[k]->stars.count))) {
+
scheduler_splittask_hydro(
scheduler_addtask(s, task_type_self, t->subtype, 0, 0,
ci->progeny[k], NULL),
s);
+ }
+ }
/* Make a task for each pair of progeny */
- for (int j = 0; j < 8; j++)
- if (ci->progeny[j] != NULL && ci->progeny[j]->hydro.count)
- for (int k = j + 1; k < 8; k++)
- if (ci->progeny[k] != NULL && ci->progeny[k]->hydro.count)
+ for (int j = 0; j < 8; j++) {
+
+ /* Do we have a non-empty progenitor? */
+ if (ci->progeny[j] != NULL &&
+ (ci->progeny[j]->hydro.count ||
+ (with_feedback && ci->progeny[j]->stars.count))) {
+
+ for (int k = j + 1; k < 8; k++) {
+
+ /* Do we have a second non-empty progenitor? */
+ if (ci->progeny[k] != NULL &&
+ (ci->progeny[k]->hydro.count ||
+ (with_feedback && ci->progeny[k]->stars.count))) {
+
scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype,
sub_sid_flag[j][k], 0, ci->progeny[j],
ci->progeny[k]),
s);
+ }
+ }
+ }
+ }
}
+
} /* Cell is split */
} /* Self interaction */
@@ -1014,478 +1052,6 @@ static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
} /* iterate over the current task. */
}
-/**
- * @brief Split a stars task if too large.
- *
- * @param t The #task
- * @param s The #scheduler we are working in.
- */
-static void scheduler_splittask_stars(struct task *t, struct scheduler *s) {
-
- /* Iterate on this task until we're done with it. */
- int redo = 1;
- while (redo) {
-
- /* Reset the redo flag. */
- redo = 0;
-
- /* Empty task? */
- if ((t->ci == NULL) || (t->type == task_type_pair && t->cj == NULL) ||
- t->ci->stars.count == 0 || (t->cj != NULL && t->cj->stars.count == 0)) {
- t->type = task_type_none;
- t->subtype = task_subtype_none;
- t->cj = NULL;
- t->skip = 1;
- break;
- }
-
- /* Self-interaction? */
- if (t->type == task_type_self) {
-
- /* Get a handle on the cell involved. */
- struct cell *ci = t->ci;
-
- /* Foreign task? */
- if (ci->nodeID != s->nodeID) {
- t->skip = 1;
- break;
- }
-
- /* Is this cell even split and the task does not violate h ? */
- if (cell_can_split_self_stars_task(ci)) {
-
- /* Make a sub? */
- if (scheduler_dosub && ci->stars.count < space_subsize_self_stars) {
-
- /* convert to a self-subtask. */
- t->type = task_type_sub_self;
-
- /* Otherwise, make tasks explicitly. */
- } else {
-
- /* Take a step back (we're going to recycle the current task)... */
- redo = 1;
-
- /* Add the self tasks. */
- int first_child = 0;
- while (ci->progeny[first_child] == NULL) first_child++;
- t->ci = ci->progeny[first_child];
- for (int k = first_child + 1; k < 8; k++)
- if (ci->progeny[k] != NULL && ci->progeny[k]->stars.count)
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_self, t->subtype, 0, 0,
- ci->progeny[k], NULL),
- s);
-
- /* Make a task for each pair of progeny */
- for (int j = 0; j < 8; j++)
- if (ci->progeny[j] != NULL && ci->progeny[j]->stars.count)
- for (int k = j + 1; k < 8; k++)
- if (ci->progeny[k] != NULL && ci->progeny[k]->stars.count)
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype,
- sub_sid_flag[j][k], 0, ci->progeny[j],
- ci->progeny[k]),
- s);
- }
- } /* Cell is split */
-
- } /* Self interaction */
-
- /* Pair interaction? */
- else if (t->type == task_type_pair) {
-
- /* Get a handle on the cells involved. */
- struct cell *ci = t->ci;
- struct cell *cj = t->cj;
-
- /* Foreign task? */
- if (ci->nodeID != s->nodeID && cj->nodeID != s->nodeID) {
- t->skip = 1;
- break;
- }
-
- /* Get the sort ID, use space_getsid and not t->flags
- to make sure we get ci and cj swapped if needed. */
- double shift[3];
- const int sid = space_getsid(s->space, &ci, &cj, shift);
-
- /* Should this task be split-up? */
- if (cell_can_split_pair_stars_task(ci) &&
- cell_can_split_pair_stars_task(cj)) {
-
- /* Replace by a single sub-task? */
- if (scheduler_dosub && /* Use division to avoid integer overflow. */
- ci->stars.count * sid_scale[sid] <
- space_subsize_pair_stars / cj->stars.count &&
- !sort_is_corner(sid)) {
-
- /* Make this task a sub task. */
- t->type = task_type_sub_pair;
-
- /* Otherwise, split it. */
- } else {
-
- /* Take a step back (we're going to recycle the current task)... */
- redo = 1;
-
- /* For each different sorting type... */
- switch (sid) {
-
- case 0: /* ( 1 , 1 , 1 ) */
- t->ci = ci->progeny[7];
- t->cj = cj->progeny[0];
- t->flags = 0;
- break;
-
- case 1: /* ( 1 , 1 , 0 ) */
- t->ci = ci->progeny[6];
- t->cj = cj->progeny[0];
- t->flags = 1;
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 1, 0,
- ci->progeny[7], cj->progeny[1]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
- ci->progeny[6], cj->progeny[1]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
- ci->progeny[7], cj->progeny[0]),
- s);
- break;
-
- case 2: /* ( 1 , 1 , -1 ) */
- t->ci = ci->progeny[6];
- t->cj = cj->progeny[1];
- t->flags = 2;
- break;
-
- case 3: /* ( 1 , 0 , 1 ) */
- t->ci = ci->progeny[5];
- t->cj = cj->progeny[0];
- t->flags = 3;
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 3, 0,
- ci->progeny[7], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
- ci->progeny[5], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 6, 0,
- ci->progeny[7], cj->progeny[0]),
- s);
- break;
-
- case 4: /* ( 1 , 0 , 0 ) */
- t->ci = ci->progeny[4];
- t->cj = cj->progeny[0];
- t->flags = 4;
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 5, 0,
- ci->progeny[5], cj->progeny[0]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
- ci->progeny[6], cj->progeny[0]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
- ci->progeny[7], cj->progeny[0]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 3, 0,
- ci->progeny[4], cj->progeny[1]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 4, 0,
- ci->progeny[5], cj->progeny[1]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 6, 0,
- ci->progeny[6], cj->progeny[1]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
- ci->progeny[7], cj->progeny[1]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 1, 0,
- ci->progeny[4], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
- ci->progeny[5], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 4, 0,
- ci->progeny[6], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 5, 0,
- ci->progeny[7], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
- ci->progeny[4], cj->progeny[3]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 1, 0,
- ci->progeny[5], cj->progeny[3]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 3, 0,
- ci->progeny[6], cj->progeny[3]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 4, 0,
- ci->progeny[7], cj->progeny[3]),
- s);
- break;
-
- case 5: /* ( 1 , 0 , -1 ) */
- t->ci = ci->progeny[4];
- t->cj = cj->progeny[1];
- t->flags = 5;
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 5, 0,
- ci->progeny[6], cj->progeny[3]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
- ci->progeny[4], cj->progeny[3]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
- ci->progeny[6], cj->progeny[1]),
- s);
- break;
-
- case 6: /* ( 1 , -1 , 1 ) */
- t->ci = ci->progeny[5];
- t->cj = cj->progeny[2];
- t->flags = 6;
- break;
-
- case 7: /* ( 1 , -1 , 0 ) */
- t->ci = ci->progeny[4];
- t->cj = cj->progeny[3];
- t->flags = 6;
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
- ci->progeny[5], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
- ci->progeny[4], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
- ci->progeny[5], cj->progeny[3]),
- s);
- break;
-
- case 8: /* ( 1 , -1 , -1 ) */
- t->ci = ci->progeny[4];
- t->cj = cj->progeny[3];
- t->flags = 8;
- break;
-
- case 9: /* ( 0 , 1 , 1 ) */
- t->ci = ci->progeny[3];
- t->cj = cj->progeny[0];
- t->flags = 9;
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 9, 0,
- ci->progeny[7], cj->progeny[4]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
- ci->progeny[3], cj->progeny[4]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
- ci->progeny[7], cj->progeny[0]),
- s);
- break;
-
- case 10: /* ( 0 , 1 , 0 ) */
- t->ci = ci->progeny[2];
- t->cj = cj->progeny[0];
- t->flags = 10;
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 11, 0,
- ci->progeny[3], cj->progeny[0]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
- ci->progeny[6], cj->progeny[0]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 6, 0,
- ci->progeny[7], cj->progeny[0]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 9, 0,
- ci->progeny[2], cj->progeny[1]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 10, 0,
- ci->progeny[3], cj->progeny[1]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
- ci->progeny[6], cj->progeny[1]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
- ci->progeny[7], cj->progeny[1]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 1, 0,
- ci->progeny[2], cj->progeny[4]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
- ci->progeny[3], cj->progeny[4]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 10, 0,
- ci->progeny[6], cj->progeny[4]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 11, 0,
- ci->progeny[7], cj->progeny[4]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
- ci->progeny[2], cj->progeny[5]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 1, 0,
- ci->progeny[3], cj->progeny[5]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 9, 0,
- ci->progeny[6], cj->progeny[5]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 10, 0,
- ci->progeny[7], cj->progeny[5]),
- s);
- break;
-
- case 11: /* ( 0 , 1 , -1 ) */
- t->ci = ci->progeny[2];
- t->cj = cj->progeny[1];
- t->flags = 11;
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 11, 0,
- ci->progeny[6], cj->progeny[5]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
- ci->progeny[2], cj->progeny[5]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 6, 0,
- ci->progeny[6], cj->progeny[1]),
- s);
- break;
-
- case 12: /* ( 0 , 0 , 1 ) */
- t->ci = ci->progeny[1];
- t->cj = cj->progeny[0];
- t->flags = 12;
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 11, 0,
- ci->progeny[3], cj->progeny[0]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 5, 0,
- ci->progeny[5], cj->progeny[0]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
- ci->progeny[7], cj->progeny[0]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 9, 0,
- ci->progeny[1], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 12, 0,
- ci->progeny[3], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
- ci->progeny[5], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 5, 0,
- ci->progeny[7], cj->progeny[2]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 3, 0,
- ci->progeny[1], cj->progeny[4]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 6, 0,
- ci->progeny[3], cj->progeny[4]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 12, 0,
- ci->progeny[5], cj->progeny[4]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 11, 0,
- ci->progeny[7], cj->progeny[4]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
- ci->progeny[1], cj->progeny[6]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 3, 0,
- ci->progeny[3], cj->progeny[6]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 9, 0,
- ci->progeny[5], cj->progeny[6]),
- s);
- scheduler_splittask_stars(
- scheduler_addtask(s, task_type_pair, t->subtype, 12, 0,
- ci->progeny[7], cj->progeny[6]),
- s);
- break;
- } /* switch(sid) */
- }
-
- /* Otherwise, break it up if it is too large? */
- } else if (scheduler_doforcesplit && ci->split && cj->split &&
- (ci->stars.count > space_maxsize / cj->stars.count)) {
-
- /* Replace the current task. */
- t->type = task_type_none;
-
- for (int j = 0; j < 8; j++)
- if (ci->progeny[j] != NULL && ci->progeny[j]->stars.count)
- for (int k = 0; k < 8; k++)
- if (cj->progeny[k] != NULL && cj->progeny[k]->stars.count) {
- struct task *tl =
- scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
- ci->progeny[j], cj->progeny[k]);
- scheduler_splittask_stars(tl, s);
- tl->flags = space_getsid(s->space, &t->ci, &t->cj, shift);
- }
- }
- } /* pair interaction? */
- } /* iterate over the current task. */
-}
-
/**
* @brief Split a gravity task if too large.
*
@@ -1670,8 +1236,6 @@ void scheduler_splittasks_mapper(void *map_data, int num_elements,
scheduler_splittask_gravity(t, s);
} else if (t->type == task_type_grav_mesh) {
/* For future use */
- } else if (t->subtype == task_subtype_stars_density) {
- scheduler_splittask_stars(t, s);
} else {
#ifdef SWIFT_DEBUG_CHECKS
error("Unexpected task sub-type");
@@ -2071,8 +1635,11 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
case task_type_grav_mm:
cost = wscale * (gcount_i + gcount_j);
break;
- case task_type_end_force:
- cost = wscale * count_i + wscale * gcount_i;
+ case task_type_end_hydro_force:
+ cost = wscale * count_i;
+ break;
+ case task_type_end_grav_force:
+ cost = wscale * gcount_i;
break;
case task_type_kick1:
cost = wscale * count_i + wscale * gcount_i;
@@ -2286,9 +1853,6 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
err = MPI_Irecv(t->ci->hydro.parts, t->ci->hydro.count, part_mpi_type,
t->ci->nodeID, t->flags, subtaskMPI_comms[t->subtype],
&t->req);
- // message( "receiving %i parts with tag=%i from %i to %i." ,
- // t->ci->hydro.count , t->flags , t->ci->nodeID , s->nodeID );
- // fflush(stdout);
} else if (t->subtype == task_subtype_gpart) {
err = MPI_Irecv(t->ci->grav.parts, t->ci->grav.count, gpart_mpi_type,
t->ci->nodeID, t->flags, subtaskMPI_comms[t->subtype],
@@ -2342,9 +1906,6 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
err = MPI_Issend(t->ci->hydro.parts, t->ci->hydro.count,
part_mpi_type, t->cj->nodeID, t->flags,
subtaskMPI_comms[t->subtype], &t->req);
- // message( "sending %i parts with tag=%i from %i to %i." ,
- // t->ci->hydro.count , t->flags , s->nodeID , t->cj->nodeID );
- // fflush(stdout);
} else if (t->subtype == task_subtype_gpart) {
if ((t->ci->grav.count * sizeof(struct gpart)) > s->mpi_message_limit)
err = MPI_Isend(t->ci->grav.parts, t->ci->grav.count,
diff --git a/src/scheduler.h b/src/scheduler.h
index f1e130c6ce2a8538b0126e86ee0cbd79cf5a0e0d..3fb1c2e43bddd73a08edd40c9c2969ea20ee6d39 100644
--- a/src/scheduler.h
+++ b/src/scheduler.h
@@ -120,6 +120,7 @@ struct scheduler {
*/
__attribute__((always_inline)) INLINE static void scheduler_activate(
struct scheduler *s, struct task *t) {
+
if (atomic_cas(&t->skip, 1, 0)) {
t->wait = 0;
int ind = atomic_inc(&s->active_count);
@@ -143,7 +144,9 @@ scheduler_activate_send(struct scheduler *s, struct link *link, int nodeID) {
struct link *l = NULL;
for (l = link; l != NULL && l->t->cj->nodeID != nodeID; l = l->next)
;
- if (l == NULL) error("Missing link to send task.");
+ if (l == NULL) {
+ error("Missing link to send task.");
+ }
scheduler_activate(s, l->t);
return l;
}
diff --git a/src/space.c b/src/space.c
index 897554225e745819f9f9b9b1fc9c1e35d988713a..0d7621d4809662ec513c2717722d13871540470e 100644
--- a/src/space.c
+++ b/src/space.c
@@ -67,8 +67,6 @@ int space_subsize_pair_hydro = space_subsize_pair_hydro_default;
int space_subsize_self_hydro = space_subsize_self_hydro_default;
int space_subsize_pair_grav = space_subsize_pair_grav_default;
int space_subsize_self_grav = space_subsize_self_grav_default;
-int space_subsize_pair_stars = space_subsize_pair_stars_default;
-int space_subsize_self_stars = space_subsize_self_stars_default;
int space_subdepth_diff_grav = space_subdepth_diff_grav_default;
int space_maxsize = space_maxsize_default;
@@ -216,8 +214,6 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->hydro.ghost_in = NULL;
c->hydro.ghost_out = NULL;
c->hydro.ghost = NULL;
- c->stars.ghost_in = NULL;
- c->stars.ghost_out = NULL;
c->stars.ghost = NULL;
c->stars.density = NULL;
c->stars.feedback = NULL;
@@ -225,8 +221,11 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->kick2 = NULL;
c->timestep = NULL;
c->timestep_limiter = NULL;
- c->end_force = NULL;
+ c->hydro.end_force = NULL;
c->hydro.drift = NULL;
+ c->stars.drift = NULL;
+ c->stars.stars_in = NULL;
+ c->stars.stars_out = NULL;
c->grav.drift = NULL;
c->grav.drift_out = NULL;
c->hydro.cooling = NULL;
@@ -234,6 +233,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->grav.down_in = NULL;
c->grav.down = NULL;
c->grav.mesh = NULL;
+ c->grav.end_force = NULL;
c->super = c;
c->hydro.super = c;
c->grav.super = c;
@@ -243,8 +243,9 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->stars.parts = NULL;
c->hydro.do_sub_sort = 0;
c->stars.do_sub_sort = 0;
- c->grav.do_sub_drift = 0;
c->hydro.do_sub_drift = 0;
+ c->grav.do_sub_drift = 0;
+ c->stars.do_sub_drift = 0;
c->hydro.do_sub_limiter = 0;
c->hydro.do_limiter = 0;
c->hydro.ti_end_min = -1;
@@ -252,6 +253,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->grav.ti_end_min = -1;
c->grav.ti_end_max = -1;
c->stars.ti_end_min = -1;
+ c->stars.ti_end_max = -1;
#ifdef SWIFT_DEBUG_CHECKS
c->cellID = 0;
#endif
@@ -274,6 +276,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->mpi.hydro.recv_rho = NULL;
c->mpi.hydro.recv_gradient = NULL;
c->mpi.grav.recv = NULL;
+ c->mpi.stars.recv = NULL;
c->mpi.recv_ti = NULL;
c->mpi.limiter.recv = NULL;
@@ -281,6 +284,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->mpi.hydro.send_rho = NULL;
c->mpi.hydro.send_gradient = NULL;
c->mpi.grav.send = NULL;
+ c->mpi.stars.send = NULL;
c->mpi.send_ti = NULL;
c->mpi.limiter.send = NULL;
#endif
@@ -539,6 +543,7 @@ void space_regrid(struct space *s, int verbose) {
c->grav.super = c;
c->hydro.ti_old_part = ti_current;
c->grav.ti_old_part = ti_current;
+ c->stars.ti_old_part = ti_current;
c->grav.ti_old_multipole = ti_current;
#ifdef WITH_MPI
c->mpi.tag = -1;
@@ -548,6 +553,8 @@ void space_regrid(struct space *s, int verbose) {
c->mpi.hydro.send_xv = NULL;
c->mpi.hydro.send_rho = NULL;
c->mpi.hydro.send_gradient = NULL;
+ c->mpi.stars.send = NULL;
+ c->mpi.stars.recv = NULL;
c->mpi.grav.recv = NULL;
c->mpi.grav.send = NULL;
#endif // WITH_MPI
@@ -1526,6 +1533,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
c->hydro.ti_old_part = ti_current;
c->grav.ti_old_part = ti_current;
c->grav.ti_old_multipole = ti_current;
+ c->stars.ti_old_part = ti_current;
#ifdef SWIFT_DEBUG_CHECKS
c->cellID = -last_cell_id;
@@ -2592,7 +2600,8 @@ void space_split_recursive(struct space *s, struct cell *c,
ti_hydro_beg_max = 0;
integertime_t ti_gravity_end_min = max_nr_timesteps, ti_gravity_end_max = 0,
ti_gravity_beg_max = 0;
- integertime_t ti_stars_end_min = max_nr_timesteps;
+ integertime_t ti_stars_end_min = max_nr_timesteps, ti_stars_end_max = 0,
+ ti_stars_beg_max = 0;
struct part *parts = c->hydro.parts;
struct gpart *gparts = c->grav.parts;
struct spart *sparts = c->stars.parts;
@@ -2685,6 +2694,7 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->hydro.ti_old_part = c->hydro.ti_old_part;
cp->grav.ti_old_part = c->grav.ti_old_part;
cp->grav.ti_old_multipole = c->grav.ti_old_multipole;
+ cp->stars.ti_old_part = c->stars.ti_old_part;
cp->loc[0] = c->loc[0];
cp->loc[1] = c->loc[1];
cp->loc[2] = c->loc[2];
@@ -2712,6 +2722,7 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->stars.do_sub_sort = 0;
cp->grav.do_sub_drift = 0;
cp->hydro.do_sub_drift = 0;
+ cp->stars.do_sub_drift = 0;
cp->hydro.do_sub_limiter = 0;
cp->hydro.do_limiter = 0;
#ifdef WITH_MPI
@@ -2762,6 +2773,8 @@ void space_split_recursive(struct space *s, struct cell *c,
ti_gravity_end_max = max(ti_gravity_end_max, cp->grav.ti_end_max);
ti_gravity_beg_max = max(ti_gravity_beg_max, cp->grav.ti_beg_max);
ti_stars_end_min = min(ti_stars_end_min, cp->stars.ti_end_min);
+ ti_stars_end_max = min(ti_stars_end_max, cp->stars.ti_end_max);
+ ti_stars_beg_max = min(ti_stars_beg_max, cp->stars.ti_beg_max);
/* Increase the depth */
if (cp->maxdepth > maxdepth) maxdepth = cp->maxdepth;
@@ -2990,6 +3003,8 @@ void space_split_recursive(struct space *s, struct cell *c,
c->grav.ti_end_max = ti_gravity_end_max;
c->grav.ti_beg_max = ti_gravity_beg_max;
c->stars.ti_end_min = ti_stars_end_min;
+ c->stars.ti_end_max = ti_stars_end_max;
+ c->stars.ti_beg_max = ti_stars_beg_max;
c->stars.h_max = stars_h_max;
c->maxdepth = maxdepth;
@@ -3509,6 +3524,8 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count,
const ptrdiff_t delta = sp - s->sparts;
#endif
+ const float initial_h = s->initial_spart_h;
+
const int with_feedback = (e->policy & engine_policy_feedback);
const struct cosmology *cosmo = e->cosmology;
@@ -3521,6 +3538,11 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count,
sp[k].v[1] *= a_factor_vel;
sp[k].v[2] *= a_factor_vel;
+ /* Imposed smoothing length from parameter file */
+ if (initial_h != -1.f) {
+ sp[k].h = initial_h;
+ }
+
#ifdef HYDRO_DIMENSION_2D
sp[k].x[2] = 0.f;
sp[k].v[2] = 0.f;
@@ -3816,12 +3838,6 @@ void space_init(struct space *s, struct swift_params *params,
space_subsize_self_grav =
parser_get_opt_param_int(params, "Scheduler:cell_sub_size_self_grav",
space_subsize_self_grav_default);
- space_subsize_pair_stars =
- parser_get_opt_param_int(params, "Scheduler:cell_sub_size_pair_stars",
- space_subsize_pair_stars_default);
- space_subsize_self_stars =
- parser_get_opt_param_int(params, "Scheduler:cell_sub_size_self_stars",
- space_subsize_self_stars_default);
space_splitsize = parser_get_opt_param_int(
params, "Scheduler:cell_split_size", space_splitsize_default);
space_subdepth_diff_grav =
@@ -3842,8 +3858,6 @@ void space_init(struct space *s, struct swift_params *params,
space_subsize_pair_hydro, space_subsize_self_hydro);
message("sub_size_pair_grav set to %d, sub_size_self_grav set to %d",
space_subsize_pair_grav, space_subsize_self_grav);
- message("sub_size_pair_stars set to %d, sub_size_self_stars set to %d",
- space_subsize_pair_stars, space_subsize_self_stars);
}
/* Apply h scaling */
@@ -3854,6 +3868,13 @@ void space_init(struct space *s, struct swift_params *params,
for (size_t k = 0; k < Npart; k++) parts[k].h *= scaling;
}
+ /* Read in imposed star smoothing length */
+ s->initial_spart_h = parser_get_opt_param_float(
+ params, "InitialConditions:stars_smoothing_length", -1.f);
+ if (s->initial_spart_h != -1.f) {
+ message("Imposing a star smoothing length of %e", s->initial_spart_h);
+ }
+
/* Apply shift */
double shift[3] = {0.0, 0.0, 0.0};
parser_get_opt_param_double_array(params, "InitialConditions:shift", 3,
@@ -4283,6 +4304,7 @@ void space_check_drift_point(struct space *s, integertime_t ti_drift,
/* Recursively check all cells */
space_map_cells_pre(s, 1, cell_check_part_drift_point, &ti_drift);
space_map_cells_pre(s, 1, cell_check_gpart_drift_point, &ti_drift);
+ space_map_cells_pre(s, 1, cell_check_spart_drift_point, &ti_drift);
if (multipole)
space_map_cells_pre(s, 1, cell_check_multipole_drift_point, &ti_drift);
#else
diff --git a/src/space.h b/src/space.h
index 98ab2523668c9789bb644f0ebe300cf73ef6f182..fe47a2b8b995e5872e79e755b5b8075a409795b8 100644
--- a/src/space.h
+++ b/src/space.h
@@ -53,8 +53,6 @@ struct cosmology;
#define space_subsize_self_hydro_default 32000
#define space_subsize_pair_grav_default 256000000
#define space_subsize_self_grav_default 32000
-#define space_subsize_pair_stars_default 256000000
-#define space_subsize_self_stars_default 32000
#define space_subdepth_diff_grav_default 4
#define space_max_top_level_cells_default 12
#define space_stretch 1.10f
@@ -226,6 +224,9 @@ struct space {
/*! Sum of the norm of the velocity of all the #spart */
float sum_spart_vel_norm;
+ /*! Initial value of the smoothing length read from the parameter file */
+ float initial_spart_h;
+
/*! General-purpose lock for this space. */
swift_lock_type lock;
diff --git a/src/stars/Default/stars.h b/src/stars/Default/stars.h
index ab4c6c7013d47a5731440bc953ad6a1101c7d2a4..586a87f75600a08acfd84b0f7ecc57fc4573281f 100644
--- a/src/stars/Default/stars.h
+++ b/src/stars/Default/stars.h
@@ -65,6 +65,26 @@ __attribute__((always_inline)) INLINE static void stars_init_spart(
sp->density.wcount_dh = 0.f;
}
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * @param sp The particle
+ * @param dt_drift The drift time-step for positions.
+ */
+__attribute__((always_inline)) INLINE static void stars_predict_extra(
+ struct spart* restrict sp, float dt_drift) {
+
+ // MATTHIEU
+ /* const float h_inv = 1.f / sp->h; */
+
+ /* /\* Predict smoothing length *\/ */
+ /* const float w1 = sp->feedback.h_dt * h_inv * dt_drift; */
+ /* if (fabsf(w1) < 0.2f) */
+ /* sp->h *= approx_expf(w1); /\* 4th order expansion of exp(w) *\/ */
+ /* else */
+ /* sp->h *= expf(w1); */
+}
+
/**
* @brief Sets the values to be predicted in the drifts to their values at a
* kick time
@@ -77,12 +97,13 @@ __attribute__((always_inline)) INLINE static void stars_reset_predicted_values(
/**
* @brief Finishes the calculation of (non-gravity) forces acting on stars
*
- * Multiplies the forces and accelerations by the appropiate constants
- *
* @param sp The particle to act upon
*/
-__attribute__((always_inline)) INLINE static void stars_end_force(
- struct spart* sp) {}
+__attribute__((always_inline)) INLINE static void stars_end_feedback(
+ struct spart* sp) {
+
+ sp->feedback.h_dt *= sp->h * hydro_dimension_inv;
+}
/**
* @brief Kick the additional variables
@@ -147,4 +168,25 @@ __attribute__((always_inline)) INLINE static void stars_evolve_spart(
struct spart* restrict sp, const struct stars_props* stars_properties,
const struct cosmology* cosmo) {}
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * This is the equivalent of hydro_reset_acceleration.
+ * We do not compute the acceleration on star, therefore no need to use it.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void stars_reset_feedback(
+ struct spart* restrict p) {
+
+ /* Reset time derivative */
+ p->feedback.h_dt = 0.f;
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+ p->ids_ngbs_force[i] = -1;
+ p->num_ngb_force = 0;
+#endif
+}
+
#endif /* SWIFT_DEFAULT_STARS_H */
diff --git a/src/stars/Default/stars_iact.h b/src/stars/Default/stars_iact.h
index 9e27f86028245a230cfd777dfc46da7b7d2f3915..994ff98d01920124081598d7e6758c929c6d0b2f 100644
--- a/src/stars/Default/stars_iact.h
+++ b/src/stars/Default/stars_iact.h
@@ -35,6 +35,8 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
/* Update ngb counters */
if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS)
si->ids_ngbs_density[si->num_ngb_density] = pj->id;
+
+ /* Update ngb counters */
++si->num_ngb_density;
#endif
}
@@ -54,4 +56,34 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
struct spart *restrict si,
- struct part *restrict pj, float a, float H) {}
+ struct part *restrict pj, float a, float H) {
+
+ const float mj = pj->mass;
+ const float rhoj = pj->rho;
+ const float r = sqrtf(r2);
+ const float ri = 1.f / r;
+
+ /* Get the kernel for hi. */
+ float hi_inv = 1.0f / hi;
+ float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+ float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+ float wi_dr = hid_inv * wi_dx;
+
+ /* Compute dv dot r */
+ float dvdr = (si->v[0] - pj->v[0]) * dx[0] + (si->v[1] - pj->v[1]) * dx[1] +
+ (si->v[2] - pj->v[2]) * dx[2];
+
+ /* Get the time derivative for h. */
+ si->feedback.h_dt -= mj * dvdr * ri / rhoj * wi_dr;
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ /* Update ngb counters */
+ if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_STARS)
+ si->ids_ngbs_force[si->num_ngb_force] = pj->id;
+
+ /* Update ngb counters */
+ ++si->num_ngb_force;
+#endif
+}
diff --git a/src/stars/Default/stars_io.h b/src/stars/Default/stars_io.h
index a6c2768f715e3dc6e870ee92e7d8a5e9458a5d11..26a42b8c0f80beb32695e2cb00716f283289663d 100644
--- a/src/stars/Default/stars_io.h
+++ b/src/stars/Default/stars_io.h
@@ -60,10 +60,10 @@ INLINE static void stars_write_particles(const struct spart *sparts,
struct io_props *list,
int *num_fields) {
- /* Say how much we want to read */
+ /* Say how much we want to write */
*num_fields = 5;
- /* List what we want to read */
+ /* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
sparts, x);
list[1] =
@@ -74,6 +74,23 @@ INLINE static void stars_write_particles(const struct spart *sparts,
sparts, id);
list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
sparts, h);
+
+#ifdef DEBUG_INTERACTIONS_STARS
+
+ list += *num_fields;
+ *num_fields += 4;
+
+ list[0] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS,
+ sparts, num_ngb_density);
+ list[1] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS,
+ sparts, num_ngb_force);
+ list[2] = io_make_output_field("Ids_ngb_density", LONGLONG,
+ MAX_NUM_OF_NEIGHBOURS_STARS,
+ UNIT_CONV_NO_UNITS, sparts, ids_ngbs_density);
+ list[3] = io_make_output_field("Ids_ngb_force", LONGLONG,
+ MAX_NUM_OF_NEIGHBOURS_STARS,
+ UNIT_CONV_NO_UNITS, sparts, ids_ngbs_force);
+#endif
}
/**
@@ -108,9 +125,6 @@ INLINE static void stars_props_init(struct stars_props *sp,
(pow_dimension(delta_eta) - pow_dimension(sp->eta_neighbours)) *
kernel_norm;
- /* Maximal smoothing length */
- sp->h_max = parser_get_opt_param_float(params, "Stars:h_max", p->h_max);
-
/* Number of iterations to converge h */
sp->max_smoothing_iterations = parser_get_opt_param_int(
params, "Stars:max_ghost_iterations", p->max_smoothing_iterations);
@@ -143,9 +157,6 @@ INLINE static void stars_props_print(const struct stars_props *sp) {
"(max|dlog(h)/dt|=%f).",
pow_dimension(expf(sp->log_max_h_change)), sp->log_max_h_change);
- if (sp->h_max != FLT_MAX)
- message("Maximal smoothing length allowed: %.4f", sp->h_max);
-
message("Maximal iterations in ghost task set to %d",
sp->max_smoothing_iterations);
}
@@ -161,7 +172,6 @@ INLINE static void stars_props_print_snapshot(hid_t h_grpstars,
io_write_attribute_f(h_grpstars, "Kernel eta", sp->eta_neighbours);
io_write_attribute_f(h_grpstars, "Smoothing length tolerance",
sp->h_tolerance);
- io_write_attribute_f(h_grpstars, "Maximal smoothing length", sp->h_max);
io_write_attribute_f(h_grpstars, "Volume log(max(delta h))",
sp->log_max_h_change);
io_write_attribute_f(h_grpstars, "Volume max change time-step",
diff --git a/src/stars/Default/stars_part.h b/src/stars/Default/stars_part.h
index eb253fd3c9f0a44efeafbfdd2b002f4b2694c4e3..bed2e14756ff2b2b83dbd1f5de821aae4ca7be51 100644
--- a/src/stars/Default/stars_part.h
+++ b/src/stars/Default/stars_part.h
@@ -61,6 +61,7 @@ struct spart {
timebin_t time_bin;
struct {
+
/* Number of neighbours. */
float wcount;
@@ -69,6 +70,13 @@ struct spart {
} density;
+ struct {
+
+ /* Change in smoothing length over time. */
+ float h_dt;
+
+ } feedback;
+
/*! Tracer structure */
struct tracers_xpart_data tracers_data;
@@ -86,11 +94,17 @@ struct spart {
#endif
#ifdef DEBUG_INTERACTIONS_STARS
+ /*! Number of interactions in the density SELF and PAIR */
+ int num_ngb_density;
+
/*! List of interacting particles in the density SELF and PAIR */
long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS_STARS];
- /*! Number of interactions in the density SELF and PAIR */
- int num_ngb_density;
+ /*! Number of interactions in the force SELF and PAIR */
+ int num_ngb_force;
+
+ /*! List of interacting particles in the force SELF and PAIR */
+ long long ids_ngbs_force[MAX_NUM_OF_NEIGHBOURS_STARS];
#endif
} SWIFT_STRUCT_ALIGN;
@@ -112,9 +126,6 @@ struct stars_props {
/*! Tolerance on neighbour number (for info only)*/
float delta_neighbours;
- /*! Maximal smoothing length */
- float h_max;
-
/*! Maximal number of iterations to converge h */
int max_smoothing_iterations;
diff --git a/src/stars/EAGLE/stars.h b/src/stars/EAGLE/stars.h
index 4acfef360cedf8993702d27a748a364288052410..ea63dd84453d7c02efc1232a2d96ef14af840b29 100644
--- a/src/stars/EAGLE/stars.h
+++ b/src/stars/EAGLE/stars.h
@@ -33,6 +33,25 @@ __attribute__((always_inline)) INLINE static float stars_compute_timestep(
return FLT_MAX;
}
+/**
+ * @brief Prepares a s-particle for its interactions
+ *
+ * @param sp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void stars_init_spart(
+ struct spart* sp) {
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+ sp->ids_ngbs_density[i] = -1;
+ sp->num_ngb_density = 0;
+#endif
+
+ sp->density.wcount = 0.f;
+ sp->density.wcount_dh = 0.f;
+ sp->rho_gas = 0.f;
+}
+
/**
* @brief Initialises the s-particles for the first time
*
@@ -46,24 +65,29 @@ __attribute__((always_inline)) INLINE static void stars_first_init_spart(
sp->time_bin = 0;
sp->birth_density = -1.f;
+ sp->birth_time = -1.f;
+
+ stars_init_spart(sp);
}
/**
- * @brief Prepares a s-particle for its interactions
+ * @brief Predict additional particle fields forward in time when drifting
*
- * @param sp The particle to act upon
+ * @param sp The particle
+ * @param dt_drift The drift time-step for positions.
*/
-__attribute__((always_inline)) INLINE static void stars_init_spart(
- struct spart* sp) {
-
-#ifdef DEBUG_INTERACTIONS_STARS
- for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
- sp->ids_ngbs_density[i] = -1;
- sp->num_ngb_density = 0;
-#endif
-
- sp->density.wcount = 0.f;
- sp->density.wcount_dh = 0.f;
+__attribute__((always_inline)) INLINE static void stars_predict_extra(
+ struct spart* restrict sp, float dt_drift) {
+
+ // MATTHIEU
+ /* const float h_inv = 1.f / sp->h; */
+
+ /* /\* Predict smoothing length *\/ */
+ /* const float w1 = sp->feedback.h_dt * h_inv * dt_drift; */
+ /* if (fabsf(w1) < 0.2f) */
+ /* sp->h *= approx_expf(w1); /\* 4th order expansion of exp(w) *\/ */
+ /* else */
+ /* sp->h *= expf(w1); */
}
/**
@@ -82,8 +106,11 @@ __attribute__((always_inline)) INLINE static void stars_reset_predicted_values(
*
* @param sp The particle to act upon
*/
-__attribute__((always_inline)) INLINE static void stars_end_force(
- struct spart* sp) {}
+__attribute__((always_inline)) INLINE static void stars_end_feedback(
+ struct spart* sp) {
+
+ sp->feedback.h_dt *= sp->h * hydro_dimension_inv;
+}
/**
* @brief Kick the additional variables
@@ -110,6 +137,7 @@ __attribute__((always_inline)) INLINE static void stars_end_density(
const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
/* Finish the calculation by inserting the missing h-factors */
+ sp->rho_gas *= h_inv_dim;
sp->density.wcount *= h_inv_dim;
sp->density.wcount_dh *= h_inv_dim_plus_one;
}
@@ -124,14 +152,10 @@ __attribute__((always_inline)) INLINE static void stars_end_density(
__attribute__((always_inline)) INLINE static void stars_spart_has_no_neighbours(
struct spart* restrict sp, const struct cosmology* cosmo) {
- /* Some smoothing length multiples. */
- const float h = sp->h;
- const float h_inv = 1.0f / h; /* 1/h */
- const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
-
/* Re-set problematic values */
- sp->density.wcount = kernel_root * h_inv_dim;
+ sp->density.wcount = 0.f;
sp->density.wcount_dh = 0.f;
+ sp->rho_gas = 0.f;
}
/**
@@ -148,4 +172,25 @@ __attribute__((always_inline)) INLINE static void stars_evolve_spart(
struct spart* restrict sp, const struct stars_props* stars_properties,
const struct cosmology* cosmo) {}
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * This is the equivalent of hydro_reset_acceleration.
+ * We do not compute the acceleration on star, therefore no need to use it.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void stars_reset_feedback(
+ struct spart* restrict p) {
+
+ /* Reset time derivative */
+ p->feedback.h_dt = 0.f;
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+ p->ids_ngbs_force[i] = -1;
+ p->num_ngb_force = 0;
+#endif
+}
+
#endif /* SWIFT_EAGLE_STARS_H */
diff --git a/src/stars/EAGLE/stars_iact.h b/src/stars/EAGLE/stars_iact.h
index 9e27f86028245a230cfd777dfc46da7b7d2f3915..5c4d94c2060ec42633a5eec472d844ca919b56eb 100644
--- a/src/stars/EAGLE/stars_iact.h
+++ b/src/stars/EAGLE/stars_iact.h
@@ -16,6 +16,9 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
const struct part *restrict pj, float a,
float H) {
+ /* Get the gas mass. */
+ const float mj = pj->mass;
+
float wi, wi_dx;
/* Get r and 1/r. */
@@ -31,6 +34,9 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
si->density.wcount += wi;
si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+ /* Compute contribution to the density */
+ si->rho_gas += mj * wi;
+
#ifdef DEBUG_INTERACTIONS_STARS
/* Update ngb counters */
if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS)
diff --git a/src/stars/EAGLE/stars_io.h b/src/stars/EAGLE/stars_io.h
index 64fde1a779966cc072a750fab89c74c5329bc03b..d93b4bf7cf81c56bb65d7d5d8801f843a492ead0 100644
--- a/src/stars/EAGLE/stars_io.h
+++ b/src/stars/EAGLE/stars_io.h
@@ -62,7 +62,7 @@ INLINE static void stars_write_particles(const struct spart *sparts,
int *num_fields) {
/* Say how much we want to write */
- *num_fields = 8;
+ *num_fields = 9;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -77,10 +77,12 @@ INLINE static void stars_write_particles(const struct spart *sparts,
sparts, h);
list[5] = io_make_output_field("BirthDensity", FLOAT, 1, UNIT_CONV_DENSITY,
sparts, birth_density);
- list[6] = io_make_output_field("Initial_Masses", FLOAT, 1, UNIT_CONV_MASS,
+ list[6] = io_make_output_field("InitialMasses", FLOAT, 1, UNIT_CONV_MASS,
sparts, mass_init);
- list[7] = io_make_output_field("Birth_time", FLOAT, 1, UNIT_CONV_TIME, sparts,
+ list[7] = io_make_output_field("BirthTime", FLOAT, 1, UNIT_CONV_TIME, sparts,
birth_time);
+ list[8] = io_make_output_field("GasDensity", FLOAT, 1, UNIT_CONV_DENSITY,
+ sparts, rho_gas);
}
/**
@@ -115,9 +117,6 @@ INLINE static void stars_props_init(struct stars_props *sp,
(pow_dimension(delta_eta) - pow_dimension(sp->eta_neighbours)) *
kernel_norm;
- /* Maximal smoothing length */
- sp->h_max = parser_get_opt_param_float(params, "Stars:h_max", p->h_max);
-
/* Number of iterations to converge h */
sp->max_smoothing_iterations = parser_get_opt_param_int(
params, "Stars:max_ghost_iterations", p->max_smoothing_iterations);
@@ -153,9 +152,6 @@ INLINE static void stars_props_print(const struct stars_props *sp) {
"(max|dlog(h)/dt|=%f).",
pow_dimension(expf(sp->log_max_h_change)), sp->log_max_h_change);
- if (sp->h_max != FLT_MAX)
- message("Maximal smoothing length allowed: %.4f", sp->h_max);
-
message("Maximal iterations in ghost task set to %d",
sp->max_smoothing_iterations);
}
@@ -171,7 +167,6 @@ INLINE static void stars_props_print_snapshot(hid_t h_grpstars,
io_write_attribute_f(h_grpstars, "Kernel eta", sp->eta_neighbours);
io_write_attribute_f(h_grpstars, "Smoothing length tolerance",
sp->h_tolerance);
- io_write_attribute_f(h_grpstars, "Maximal smoothing length", sp->h_max);
io_write_attribute_f(h_grpstars, "Volume log(max(delta h))",
sp->log_max_h_change);
io_write_attribute_f(h_grpstars, "Volume max change time-step",
diff --git a/src/stars/EAGLE/stars_part.h b/src/stars/EAGLE/stars_part.h
index a18a55e959078201ada3c7589db92e5b4946ad25..664b5c0d03dd5b762aab0a0f582f22312b2196c6 100644
--- a/src/stars/EAGLE/stars_part.h
+++ b/src/stars/EAGLE/stars_part.h
@@ -58,13 +58,17 @@ struct spart {
/*! Initial star mass */
float mass_init;
- /* Particle cutoff radius. */
+ /*! Particle smoothing length. */
float h;
+ /*! Density of the gas surrounding the star. */
+ float rho_gas;
+
/*! Particle time bin */
timebin_t time_bin;
struct {
+
/* Number of neighbours. */
float wcount;
@@ -73,8 +77,16 @@ struct spart {
} density;
+ struct {
+
+ /* Change in smoothing length over time. */
+ float h_dt;
+
+ } feedback;
+
/*! Union for the birth time and birht scale factor */
union {
+
/*! Birth time */
float birth_time;
@@ -128,9 +140,6 @@ struct stars_props {
/*! Tolerance on neighbour number (for info only)*/
float delta_neighbours;
- /*! Maximal smoothing length */
- float h_max;
-
/*! Maximal number of iterations to converge h */
int max_smoothing_iterations;
diff --git a/src/task.c b/src/task.c
index 69837548be104f33891f5bef97f8f5c499b62812..34c636b48ed6ff3fefdf1e7847a67ca56ea79c89 100644
--- a/src/task.c
+++ b/src/task.c
@@ -61,9 +61,10 @@ const char *taskID_names[task_type_count] = {"none",
"ghost_out",
"extra_ghost",
"drift_part",
+ "drift_spart",
"drift_gpart",
"drift_gpart_out",
- "end_force",
+ "end_hydro_force",
"kick1",
"kick2",
"timestep",
@@ -75,9 +76,12 @@ const char *taskID_names[task_type_count] = {"none",
"grav_down_in",
"grav_down",
"grav_mesh",
+ "grav_end_force",
"cooling",
"star_formation",
"logger",
+ "stars_in",
+ "stars_out",
"stars_ghost_in",
"stars_ghost",
"stars_ghost_out",
@@ -144,12 +148,14 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_type_extra_ghost:
case task_type_timestep_limiter:
case task_type_cooling:
+ case task_type_end_hydro_force:
return task_action_part;
break;
case task_type_star_formation:
return task_action_all;
+ case task_type_drift_spart:
case task_type_stars_ghost:
case task_type_stars_sort:
return task_action_spart;
@@ -179,13 +185,15 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
break;
default:
- error("Unknow task_action for task");
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Unknown task_action for task %s/%s", taskID_names[t->type],
+ subtaskID_names[t->subtype]);
+#endif
return task_action_none;
break;
}
break;
- case task_type_end_force:
case task_type_kick1:
case task_type_kick2:
case task_type_logger:
@@ -198,8 +206,11 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
return task_action_part;
else if (t->ci->grav.count > 0)
return task_action_gpart;
- else
+ else {
+#ifdef SWIFT_DEBUG_CHECKS
error("Task without particles");
+#endif
+ }
break;
case task_type_init_grav:
@@ -210,18 +221,25 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_type_drift_gpart:
case task_type_grav_down:
+ case task_type_end_grav_force:
case task_type_grav_mesh:
return task_action_gpart;
break;
default:
- error("Unknown task_action for task");
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Unknown task_action for task %s/%s", taskID_names[t->type],
+ subtaskID_names[t->subtype]);
+#endif
return task_action_none;
break;
}
- /* Silence compiler warnings */
- error("Unknown task_action for task");
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Unknown task_action for task %s/%s", taskID_names[t->type],
+ subtaskID_names[t->subtype]);
+#endif
+ /* Silence compiler warnings. We should never get here. */
return task_action_none;
}
@@ -265,6 +283,8 @@ float task_overlap(const struct task *restrict ta,
if (tb->ci != NULL) size_union += tb->ci->hydro.count;
if (tb->cj != NULL) size_union += tb->cj->hydro.count;
+ if (size_union == 0) return 0.f;
+
/* Compute the intersection of the cell data. */
const size_t size_intersect = task_cell_overlap_part(ta->ci, tb->ci) +
task_cell_overlap_part(ta->ci, tb->cj) +
@@ -284,6 +304,8 @@ float task_overlap(const struct task *restrict ta,
if (tb->ci != NULL) size_union += tb->ci->grav.count;
if (tb->cj != NULL) size_union += tb->cj->grav.count;
+ if (size_union == 0) return 0.f;
+
/* Compute the intersection of the cell data. */
const size_t size_intersect = task_cell_overlap_gpart(ta->ci, tb->ci) +
task_cell_overlap_gpart(ta->ci, tb->cj) +
@@ -303,6 +325,8 @@ float task_overlap(const struct task *restrict ta,
if (tb->ci != NULL) size_union += tb->ci->stars.count;
if (tb->cj != NULL) size_union += tb->cj->stars.count;
+ if (size_union == 0) return 0.f;
+
/* Compute the intersection of the cell data. */
const size_t size_intersect = task_cell_overlap_spart(ta->ci, tb->ci) +
task_cell_overlap_spart(ta->ci, tb->cj) +
@@ -330,7 +354,6 @@ void task_unlock(struct task *t) {
/* Act based on task type. */
switch (type) {
- case task_type_end_force:
case task_type_kick1:
case task_type_kick2:
case task_type_logger:
@@ -342,12 +365,14 @@ void task_unlock(struct task *t) {
case task_type_drift_part:
case task_type_sort:
case task_type_ghost:
+ case task_type_end_hydro_force:
case task_type_timestep_limiter:
cell_unlocktree(ci);
break;
case task_type_drift_gpart:
case task_type_grav_mesh:
+ case task_type_end_grav_force:
cell_gunlocktree(ci);
break;
@@ -453,7 +478,6 @@ int task_lock(struct task *t) {
#endif
break;
- case task_type_end_force:
case task_type_kick1:
case task_type_kick2:
case task_type_logger:
@@ -469,6 +493,7 @@ int task_lock(struct task *t) {
case task_type_drift_part:
case task_type_sort:
case task_type_ghost:
+ case task_type_end_hydro_force:
case task_type_timestep_limiter:
if (ci->hydro.hold) return 0;
if (cell_locktree(ci) != 0) return 0;
@@ -480,6 +505,7 @@ int task_lock(struct task *t) {
break;
case task_type_drift_gpart:
+ case task_type_end_grav_force:
case task_type_grav_mesh:
if (ci->grav.phold) return 0;
if (cell_glocktree(ci) != 0) return 0;
@@ -655,7 +681,11 @@ void task_get_group_name(int type, int subtype, char *cluster) {
strcpy(cluster, "Density");
break;
case task_subtype_gradient:
- strcpy(cluster, "Gradient");
+ if (type == task_type_send || type == task_type_recv) {
+ strcpy(cluster, "None");
+ } else {
+ strcpy(cluster, "Gradient");
+ }
break;
case task_subtype_force:
strcpy(cluster, "Force");
@@ -667,7 +697,10 @@ void task_get_group_name(int type, int subtype, char *cluster) {
strcpy(cluster, "Timestep_limiter");
break;
case task_subtype_stars_density:
- strcpy(cluster, "Stars");
+ strcpy(cluster, "StarsDensity");
+ break;
+ case task_subtype_stars_feedback:
+ strcpy(cluster, "StarsFeedback");
break;
default:
strcpy(cluster, "None");
diff --git a/src/task.h b/src/task.h
index 35b46bd383221d767e313fa38d838777e61f0a99..704d1a5ef80f1208bce69d0acf7625fb36fa19e1 100644
--- a/src/task.h
+++ b/src/task.h
@@ -52,9 +52,10 @@ enum task_types {
task_type_ghost_out, /* Implicit */
task_type_extra_ghost,
task_type_drift_part,
+ task_type_drift_spart,
task_type_drift_gpart,
task_type_drift_gpart_out, /* Implicit */
- task_type_end_force,
+ task_type_end_hydro_force,
task_type_kick1,
task_type_kick2,
task_type_timestep,
@@ -66,12 +67,15 @@ enum task_types {
task_type_grav_down_in, /* Implicit */
task_type_grav_down,
task_type_grav_mesh,
+ task_type_end_grav_force,
task_type_cooling,
task_type_star_formation,
task_type_logger,
- task_type_stars_ghost_in,
+ task_type_stars_in, /* Implicit */
+ task_type_stars_out, /* Implicit */
+ task_type_stars_ghost_in, /* Implicit */
task_type_stars_ghost,
- task_type_stars_ghost_out,
+ task_type_stars_ghost_out, /* Implicit */
task_type_stars_sort,
task_type_count
} __attribute__((packed));
diff --git a/src/timestep.h b/src/timestep.h
index e9943a41a0536b65944f0256c827d43386aadd88..b98ce06d5e69a2e2b5cb8503322c025dc69f92c7 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -201,8 +201,14 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
new_dt_self = gravity_compute_timestep_self(
sp->gpart, a_hydro, e->gravity_properties, e->cosmology);
+ /* Limit change in smoothing length */
+ const float dt_h_change = (sp->feedback.h_dt != 0.0f)
+ ? fabsf(e->stars_properties->log_max_h_change *
+ sp->h / sp->feedback.h_dt)
+ : FLT_MAX;
+
/* Take the minimum of all */
- float new_dt = min3(new_dt_stars, new_dt_self, new_dt_ext);
+ float new_dt = min4(new_dt_stars, new_dt_self, new_dt_ext, dt_h_change);
/* Apply the maximal displacement constraint (FLT_MAX if non-cosmological)*/
new_dt = min(new_dt, e->dt_max_RMS_displacement);
@@ -212,9 +218,10 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
- if (new_dt < e->dt_min)
+ if (new_dt < e->dt_min) {
error("spart (id=%lld) wants a time-step (%e) below dt_min (%e)", sp->id,
new_dt, e->dt_min);
+ }
/* Convert to integer time */
const integertime_t new_dti = make_integer_timestep(
diff --git a/tools/plot_task_dependencies.py b/tools/plot_task_dependencies.py
index fd7c8f8f7c4165f38947828f64a7f3fd26f41ee5..14ba7c99f621c0c62c3c9cb8e8d3b0c78e491401 100644
--- a/tools/plot_task_dependencies.py
+++ b/tools/plot_task_dependencies.py
@@ -167,6 +167,8 @@ def taskIsHydro(name):
return True
if "rho" in name:
return True
+ if "gradient" in name:
+ return True
if "force" in name:
return True
if "xv" in name:
@@ -177,6 +179,9 @@ def taskIsHydro(name):
"ghost_in",
"ghost",
"ghost_out",
+ "extra_ghost",
+ "cooling",
+ "star_formation"
]
if name in task_name:
return True
diff --git a/tools/task_plots/analyse_tasks.py b/tools/task_plots/analyse_tasks.py
index e897424a95be8937073bd16adf108fa4fa1456ad..fc9df0e4797cfb16e883df551af30dc0d3244edc 100755
--- a/tools/task_plots/analyse_tasks.py
+++ b/tools/task_plots/analyse_tasks.py
@@ -76,9 +76,10 @@ TASKTYPES = [
"ghost_out",
"extra_ghost",
"drift_part",
+ "drift_spart",
"drift_gpart",
"drift_gpart_out",
- "end_force",
+ "hydro_end_force",
"kick1",
"kick2",
"timestep",
@@ -90,9 +91,12 @@ TASKTYPES = [
"grav_down_in",
"grav_down",
"grav_mesh",
+ "grav_end_force",
"cooling",
"star_formation",
"logger",
+ "stars_in",
+ "stars_out",
"stars_ghost_in",
"stars_ghost",
"stars_ghost_out",
diff --git a/tools/task_plots/plot_tasks.py b/tools/task_plots/plot_tasks.py
index 12fd4d241a268c9d45fd72f5cdda2727221ba94d..54f34b2f828895d894b84253e366173827c03158 100755
--- a/tools/task_plots/plot_tasks.py
+++ b/tools/task_plots/plot_tasks.py
@@ -161,9 +161,10 @@ TASKTYPES = [
"ghost_out",
"extra_ghost",
"drift_part",
+ "drift_spart",
"drift_gpart",
"drift_gpart_out",
- "end_force",
+ "hydro_end_force",
"kick1",
"kick2",
"timestep",
@@ -175,9 +176,12 @@ TASKTYPES = [
"grav_down_in",
"grav_down",
"grav_mesh",
+ "grav_end_force",
"cooling",
"star_formation",
"logger",
+ "stars_in",
+ "stars_out",
"stars_ghost_in",
"stars_ghost",
"stars_ghost_out",
@@ -232,6 +236,8 @@ FULLTYPES = [
"send/tend",
"recv/gpart",
"send/gpart",
+ "recv/spart",
+ "send/spart",
"self/stars_density",
"pair/stars_density",
"sub_self/stars_density",