diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst
index eab1e670cd8e3dcdef153b2ce5b3e3c58565df14..726631b7e0de1e2ed134fbba456f0a7f3288ea38 100644
--- a/doc/RTD/source/ParameterFiles/parameter_description.rst
+++ b/doc/RTD/source/ParameterFiles/parameter_description.rst
@@ -373,6 +373,15 @@ prevent the smoothing length from going below this value in dense
environments. This will lead to smoothing over more particles than specified
by :math:`\eta`.
+The optional parameter ``particle_splitting`` (Default: 0) activates the
+splitting of overly massive particles into 2. By switching this on, the code
+will loop over all the particles at every tree rebuild and split the particles
+with a mass above a fixed threshold into two copies that are slightly shifted
+(by a randomly orientated vector of norm :math:`0.2h`). Their masses and other
+relevant particle-carried quantities are then halved. The mass threshold for
+splitting is set by the parameter ``particle_splitting_mass_threshold`` which is
+specified using the internal unit system.
+
The final set of parameters in this section determine the initial and minimum
temperatures of the particles.
@@ -399,14 +408,17 @@ The full section to start a typical cosmological run would be:
.. code:: YAML
SPH:
- resolution_eta: 1.2
- CFL_condition: 0.1
- h_tolerance: 1e-4
- h_min_ratio: 0.1
- initial_temperature: 273
- minimal_temperature: 100
- H_mass_fraction: 0.755
- H_ionization_temperature: 1e4
+ resolution_eta: 1.2
+ CFL_condition: 0.1
+ h_tolerance: 1e-4
+ h_min_ratio: 0.1
+ h_max: 1. # U_L
+ initial_temperature: 273 # U_T
+ minimal_temperature: 100 # U_T
+ H_mass_fraction: 0.755
+ H_ionization_temperature: 1e4 # U_T
+ particle_splitting: 1
+ particle_splitting_mass_threshold: 5e-3 # U_M
.. _Parameters_Stars:
diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index 32399918a50782fcdac999a3efa2edce94981640..864df962f0561a3a2284c75ccc485dbb2cafc4fe 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -48,12 +48,14 @@ Gravity:
# 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).
- h_min_ratio: 0.1 # Minimal smoothing in units of softening.
- h_max: 0.5 # Maximal softening in co-moving internal units.
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- minimal_temperature: 100.0 # (internal units)
- initial_temperature: 268.7
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ h_min_ratio: 0.1 # Minimal smoothing in units of softening.
+ h_max: 0.5 # Maximal softening in co-moving internal units.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100.0 # (internal units)
+ initial_temperature: 268.7 # (internal units)
+ particle_splitting: 1 # Particle splitting is ON
+ particle_splitting_mass_threshold: 7e-4 # (internal units, i.e. 7e6 Msun ~ 4x initial gas particle mass)
# Parameters of the stars neighbour search
Stars:
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index 69b93c8e364f7260e090a7fb011f77c2d71c2f17..9c99a164b364fa719c42f7d56d2c573a6190b26d 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -48,12 +48,14 @@ Gravity:
# 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).
- h_min_ratio: 0.1 # Minimal smoothing in units of softening.
- h_max: 0.5 # Maximal softening in co-moving internal units.
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- minimal_temperature: 100.0 # (internal units)
- initial_temperature: 268.7
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ h_min_ratio: 0.1 # Minimal smoothing in units of softening.
+ h_max: 0.5 # Maximal softening in co-moving internal units.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100.0 # (internal units)
+ initial_temperature: 268.7 # (internal units)
+ particle_splitting: 1 # Particle splitting is ON
+ particle_splitting_mass_threshold: 7e-4 # (internal units, i.e. 7e6 Msun ~ 4x initial gas particle mass)
# Parameters of the stars neighbour search
Stars:
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index f9c6867139cfb835e431d50a6ae9cc3b7d242711..3d75c944396bcc6b667a830047a9eeb24cc9cdff 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -48,12 +48,14 @@ Gravity:
# 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).
- h_min_ratio: 0.1 # Minimal smoothing in units of softening.
- h_max: 0.5 # Maximal softening in co-moving internal units.
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- minimal_temperature: 100.0 # (internal units)
- initial_temperature: 268.7
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ h_min_ratio: 0.1 # Minimal smoothing in units of softening.
+ h_max: 0.5 # Maximal softening in co-moving internal units.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100.0 # (internal units)
+ initial_temperature: 268.7 # (internal units)
+ particle_splitting: 1 # Particle splitting is ON
+ particle_splitting_mass_threshold: 7e-4 # (internal units, i.e. 7e6 Msun ~ 4x initial gas particle mass)
# Parameters of the stars neighbour search
Stars:
diff --git a/examples/EAGLE_ICs/README b/examples/EAGLE_ICs/README
index 82d449668c9bc30e765cddc29b8615b9b612430a..ea352da4e1799bedfad89b9a7874890eef5ab0d2 100644
--- a/examples/EAGLE_ICs/README
+++ b/examples/EAGLE_ICs/README
@@ -11,6 +11,9 @@ the following changes have been made:
been changed to 1.79 ckpc. This follows the recommendations of
Ludlow et al. 2019. Old values were 0.7 pkpc and 2.69 ckpc for
all the particle species.
+ - SPH particles with a mass larger than 7*10^6 Msun (~4x the initial
+ gas particle mass) are now split into 2 equal mass particles
+ within the smoothing length of the original particle.
- The metallicity-dependent density threshold for star formation
uses the smoothed metallicities and not the raw metallicities
any more.
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index 3caab833bebb0f4d7baf306f2acb58a94c9d2fd0..75aa30efb8dc235124d05e971173c460f4829bf2 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -54,11 +54,13 @@ Gravity:
# 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).
- h_min_ratio: 0.1 # Minimal smoothing in units of softening.
- h_max: 0.5 # Maximal softening in co-moving internal units.
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- minimal_temperature: 100 # (internal units)
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ h_min_ratio: 0.1 # Minimal smoothing in units of softening.
+ h_max: 0.5 # Maximal softening in co-moving internal units.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100 # (internal units)
+ particle_splitting: 1
+ particle_splitting_mass_threshold: 7e-4 # Internal units (i.e. 7e6 Msun ~ 4 times the initial gas mass)
# Parameters of the stars neighbour search
Stars:
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index 2c955661478297c741dc0779190252d8b81c7929..c6c606c87ad81d08187e17deac0f88119381c3df 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -62,11 +62,13 @@ Gravity:
# 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).
- h_min_ratio: 0.1 # Minimal smoothing in units of softening.
- h_max: 0.5 # Maximal softening in co-moving internal units.
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- minimal_temperature: 100 # (internal units)
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ h_min_ratio: 0.1 # Minimal smoothing in units of softening.
+ h_max: 0.5 # Maximal softening in co-moving internal units.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100 # (internal units)
+ particle_splitting: 1
+ particle_splitting_mass_threshold: 7e-4 # Internal units (i.e. 7e6 Msun ~ 4 times the initial gas mass)
# Parameters of the stars neighbour search
Stars:
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index 49c631a3d931eaeeaf3bbb63593c69afd8522f62..2e0f7c0a9dd752738e107f089c8b01d890916b17 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -53,11 +53,13 @@ Gravity:
# 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).
- h_min_ratio: 0.1 # Minimal smoothing in units of softening.
- h_max: 0.5 # Maximal softening in co-moving internal units.
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- minimal_temperature: 100 # (internal units)
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ h_min_ratio: 0.1 # Minimal smoothing in units of softening.
+ h_max: 0.5 # Maximal softening in co-moving internal units.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100 # (internal units)
+ particle_splitting: 1
+ particle_splitting_mass_threshold: 7e-4 # Internal units (i.e. 7e6 Msun ~ 4 times the initial gas mass)
# Parameters of the stars neighbour search
Stars:
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index 7b9608b6021bee3687b3585a3a9b02c8773a14cc..835664cb6250779a9d9acfc03d3c643bec399610 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -63,11 +63,13 @@ Gravity:
# 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).
- h_min_ratio: 0.1 # Minimal smoothing in units of softening.
- h_max: 0.5 # Maximal softening in co-moving internal units.
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- minimal_temperature: 100 # (internal units)
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ h_min_ratio: 0.1 # Minimal smoothing in units of softening.
+ h_max: 0.5 # Maximal softening in co-moving internal units.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100 # (internal units)
+ particle_splitting: 1
+ particle_splitting_mass_threshold: 7e-4 # Internal units (i.e. 7e6 Msun ~ 4 times the initial gas mass)
# Parameters of the stars neighbour search
Stars:
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index fcba1f33795edcd1de2b12b123016f10f218ad7f..754b0d998607b1c1ebf747bcac4406f65d7e9c16 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -28,26 +28,28 @@ Cosmology:
# 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).
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- use_mass_weighted_num_ngb: 0 # (Optional) Are we using the mass-weighted definition of the number of neighbours?
- h_tolerance: 1e-4 # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
- h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
- h_min_ratio: 0. # (Optional) Minimal allowed smoothing length in units of the softening. Defaults to 0 if unspecified.
- max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step.
- max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
- initial_temperature: 0 # (Optional) Initial temperature (in internal units) to set the gas particles at start-up. Value is ignored if set to 0.
- minimal_temperature: 0 # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
- H_mass_fraction: 0.755 # (Optional) Hydrogen mass fraction used for initial conversion from temp to internal energy. Default value is derived from the physical constants.
- H_ionization_temperature: 1e4 # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas.
- viscosity_alpha: 0.8 # (Optional) Override for the initial value of the artificial viscosity. In schemes that have a fixed AV, this remains as alpha throughout the run.
- viscosity_alpha_max: 2.0 # (Optional) Maximal value for the artificial viscosity in schemes that allow alpha to vary.
- viscosity_alpha_min: 0.1 # (Optional) Minimal value for the artificial viscosity in schemes that allow alpha to vary.
- viscosity_length: 0.1 # (Optional) Decay length for the artificial viscosity in schemes that allow alpha to vary.
- diffusion_alpha: 0.0 # (Optional) Override the initial value for the thermal diffusion coefficient in schemes with thermal diffusion.
- diffusion_beta: 0.01 # (Optional) Override the decay/rise rate tuning parameter for the thermal diffusion.
- diffusion_alpha_max: 1.0 # (Optional) Override the maximal thermal diffusion coefficient that is allowed for a given particle.
- diffusion_alpha_min: 0.0 # (Optional) Override the minimal thermal diffusion coefficient that is allowed for a given particle.
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ use_mass_weighted_num_ngb: 0 # (Optional) Are we using the mass-weighted definition of the number of neighbours?
+ h_tolerance: 1e-4 # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
+ h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
+ h_min_ratio: 0. # (Optional) Minimal allowed smoothing length in units of the softening. Defaults to 0 if unspecified.
+ max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step.
+ max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
+ particle_splitting: 1 # (Optional) Are we splitting particles that are too massive (default: 0)
+ particle_splitting_mass_threshold: 7e-4 # (Optional) Mass threshold for particle splitting (in internal units)
+ initial_temperature: 0 # (Optional) Initial temperature (in internal units) to set the gas particles at start-up. Value is ignored if set to 0.
+ minimal_temperature: 0 # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
+ H_mass_fraction: 0.755 # (Optional) Hydrogen mass fraction used for initial conversion from temp to internal energy. Default value is derived from the physical constants.
+ H_ionization_temperature: 1e4 # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas.
+ viscosity_alpha: 0.8 # (Optional) Override for the initial value of the artificial viscosity. In schemes that have a fixed AV, this remains as alpha throughout the run.
+ viscosity_alpha_max: 2.0 # (Optional) Maximal value for the artificial viscosity in schemes that allow alpha to vary.
+ viscosity_alpha_min: 0.1 # (Optional) Minimal value for the artificial viscosity in schemes that allow alpha to vary.
+ viscosity_length: 0.1 # (Optional) Decay length for the artificial viscosity in schemes that allow alpha to vary.
+ diffusion_alpha: 0.0 # (Optional) Override the initial value for the thermal diffusion coefficient in schemes with thermal diffusion.
+ diffusion_beta: 0.01 # (Optional) Override the decay/rise rate tuning parameter for the thermal diffusion.
+ diffusion_alpha_max: 1.0 # (Optional) Override the maximal thermal diffusion coefficient that is allowed for a given particle.
+ diffusion_alpha_min: 0.0 # (Optional) Override the minimal thermal diffusion coefficient that is allowed for a given particle.
# Parameters of the stars neighbour search
Stars:
diff --git a/src/Makefile.am b/src/Makefile.am
index 5ecfac826c10568a0e98c7f4feaced9e9152aaa1..04f9e7c58ed10435f7ba3f636ca841fdd81c40ab 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -76,7 +76,7 @@ AM_SOURCES = space.c runner_main.c runner_doiact_hydro.c runner_doiact_limiter.c
runner_doiact_stars.c runner_doiact_black_holes.c runner_ghost.c runner_recv.c \
runner_sort.c runner_drift.c runner_black_holes.c runner_time_integration.c \
runner_doiact_hydro_vec.c runner_others.c\
- queue.c task.c cell.c engine.c engine_maketasks.c \
+ queue.c task.c cell.c engine.c engine_maketasks.c engine_split_particles.c \
engine_marktasks.c engine_drift.c engine_unskip.c engine_collect_end_of_step.c \
engine_redistribute.c engine_fof.c serial_io.c timers.c debug.c scheduler.c \
proxy.c parallel_io.c units.c common_io.c single_io.c multipole.c version.c map.c \
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index 7dda6e05261feef435d2416f0845c002a9c7a15d..d2f475af9053d148345350cc02004cf501b99165 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -403,6 +403,21 @@ __attribute__((always_inline)) INLINE static void chemistry_add_bpart_to_bpart(
bp_data->iron_mass_from_SNIa += swallowed_data->iron_mass_from_SNIa;
}
+/**
+ * @brief Split the metal content of a particle into n pieces
+ *
+ * We only need to split the fields that are not fractions.
+ *
+ * @param p The #part.
+ * @param n The number of pieces to split into.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_split_part(
+ struct part* p, const double n) {
+ p->chemistry_data.mass_from_SNIa /= n;
+ p->chemistry_data.mass_from_SNII /= n;
+ p->chemistry_data.mass_from_AGB /= n;
+}
+
/**
* @brief Returns the total metallicity (metal mass fraction) of the
* star particle to be used in feedback/enrichment related routines.
diff --git a/src/chemistry/GEAR/chemistry.h b/src/chemistry/GEAR/chemistry.h
index af8d70030f31cfa4c8a9d1ff92865e8e1405c45a..d301661a96e974e1edc2c5dd1838f8300e6e309e 100644
--- a/src/chemistry/GEAR/chemistry.h
+++ b/src/chemistry/GEAR/chemistry.h
@@ -229,7 +229,9 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
*/
__attribute__((always_inline)) INLINE static void chemistry_bpart_from_part(
struct chemistry_bpart_data* bp_data,
- const struct chemistry_part_data* p_data, const double gas_mass) {}
+ const struct chemistry_part_data* p_data, const double gas_mass) {
+ error("Loic: to be implemented");
+}
/**
* @brief Add the chemistry data of a gas particle to a black hole.
@@ -242,7 +244,9 @@ __attribute__((always_inline)) INLINE static void chemistry_bpart_from_part(
*/
__attribute__((always_inline)) INLINE static void chemistry_add_part_to_bpart(
struct chemistry_bpart_data* bp_data,
- const struct chemistry_part_data* p_data, const double gas_mass) {}
+ const struct chemistry_part_data* p_data, const double gas_mass) {
+ error("Loic: to be implemented");
+}
/**
* @brief Add the chemistry data of a black hole to another one.
@@ -254,6 +258,19 @@ __attribute__((always_inline)) INLINE static void chemistry_add_part_to_bpart(
*/
__attribute__((always_inline)) INLINE static void chemistry_add_bpart_to_bpart(
struct chemistry_bpart_data* bp_data,
- const struct chemistry_bpart_data* swallowed_data) {}
+ const struct chemistry_bpart_data* swallowed_data) {
+ error("Loic: to be implemented");
+}
+
+/**
+ * @brief Split the metal content of a particle into n pieces
+ *
+ * @param p The #part.
+ * @param n The number of pieces to split into.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_split_part(
+ struct part* p, const double n) {
+ error("Loic: to be implemented");
+}
#endif /* SWIFT_CHEMISTRY_GEAR_H */
diff --git a/src/chemistry/none/chemistry.h b/src/chemistry/none/chemistry.h
index 34db9d0d9b9213a6458d65009ff6f79644ac5a1d..4070c4fe7e942778ece03c3d77a9110e33add6eb 100644
--- a/src/chemistry/none/chemistry.h
+++ b/src/chemistry/none/chemistry.h
@@ -214,6 +214,17 @@ __attribute__((always_inline)) INLINE static void chemistry_add_bpart_to_bpart(
struct chemistry_bpart_data* bp_data,
const struct chemistry_bpart_data* swallowed_data) {}
+/**
+ * @brief Split the metal content of a particle into n pieces
+ *
+ * Nothing to do here.
+ *
+ * @param p The #part.
+ * @param n The number of pieces to split into.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_split_part(
+ struct part* p, const double n) {}
+
/**
* @brief Returns the total metallicity (metal mass fraction) of the
* star particle to be used in feedback/enrichment related routines.
diff --git a/src/cooling/Compton/cooling.h b/src/cooling/Compton/cooling.h
index e04fc6f9adaef6bd238fe32d21526db1dc282aaa..3e2c0ea8a3fe32e294152b725d1f0ab32d0e6e75 100644
--- a/src/cooling/Compton/cooling.h
+++ b/src/cooling/Compton/cooling.h
@@ -308,6 +308,19 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
return xp->cooling_data.radiated_energy;
}
+/**
+ * @brief Split the coolong content of a particle into n pieces
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+static INLINE void cooling_split_part(struct part* p, struct xpart* xp,
+ double n) {
+
+ xp->cooling_data.radiated_energy /= n;
+}
+
/**
* @brief Initialises the cooling properties.
*
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index a228eb6ee74e301d2ec6ac461ab9b7420fe66250..694d429db062ded86bf58809d96c1afb9808fba1 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -642,6 +642,18 @@ __attribute__((always_inline)) INLINE float cooling_get_radiated_energy(
return xp->cooling_data.radiated_energy;
}
+/**
+ * @brief Split the coolong content of a particle into n pieces
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+void cooling_split_part(struct part *p, struct xpart *xp, double n) {
+
+ xp->cooling_data.radiated_energy /= n;
+}
+
/**
* @brief Inject a fixed amount of energy to each particle in the simulation
* to mimic Hydrogen reionization.
diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h
index f154df3136558931a917aab2bcc27167d2a40508..c115a8ce3f3a7dbd438f4332711421f0fea132fb 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/src/cooling/EAGLE/cooling.h
@@ -71,6 +71,8 @@ float cooling_get_temperature(
float cooling_get_radiated_energy(const struct xpart *restrict xp);
+void cooling_split_part(struct part *p, struct xpart *xp, double n);
+
void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling,
const struct cosmology *cosmo,
struct space *s);
diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h
index 386ad3e0013a5884b12e589b3d1b1a09bfbb89dc..e78badfa856f1d53c36730100cfb837ac7ab6b56 100644
--- a/src/cooling/const_du/cooling.h
+++ b/src/cooling/const_du/cooling.h
@@ -226,6 +226,19 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
return xp->cooling_data.radiated_energy;
}
+/**
+ * @brief Split the coolong content of a particle into n pieces
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+static INLINE void cooling_split_part(struct part* p, struct xpart* xp,
+ double n) {
+
+ xp->cooling_data.radiated_energy /= n;
+}
+
/**
* @brief Initialises the cooling function properties from the parameter file
*
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index 3e3a44249d152e407bf9d3fdf294dd504cee2257..7820a3a5e44ad6017e8cffb24d4453988ebb6c4d 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -302,6 +302,19 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
return xp->cooling_data.radiated_energy;
}
+/**
+ * @brief Split the coolong content of a particle into n pieces
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+static INLINE void cooling_split_part(struct part* p, struct xpart* xp,
+ double n) {
+
+ xp->cooling_data.radiated_energy /= n;
+}
+
/**
* @brief Initialises the cooling properties.
*
diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h
index dec0e80ffa834e7749c8c1815a8a7441befc251b..2b646aedfebe346f366565ca14f91784f3781fa0 100644
--- a/src/cooling/grackle/cooling.h
+++ b/src/cooling/grackle/cooling.h
@@ -763,6 +763,19 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
return FLT_MAX;
}
+/**
+ * @brief Split the coolong content of a particle into n pieces
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+INLINE static void cooling_split_part(struct part* p, struct xpart* xp,
+ double n) {
+
+ error("Loic: to be implemented");
+}
+
/**
* @brief Initialises the cooling unit system.
*
diff --git a/src/cooling/none/cooling.h b/src/cooling/none/cooling.h
index 5630914dea5a3ce205c14ba1fb7a353d784df37f..27c7451517fc07c17486da0ba0ee0e0691929a54 100644
--- a/src/cooling/none/cooling.h
+++ b/src/cooling/none/cooling.h
@@ -176,6 +176,18 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
return 0.f;
}
+/**
+ * @brief Split the coolong content of a particle into n pieces
+ *
+ * Nothing to do here.
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+static INLINE void cooling_split_part(struct part* p, struct xpart* xp,
+ double n) {}
+
/**
* @brief Initialises the cooling properties.
*
diff --git a/src/engine.c b/src/engine.c
index 737a84daecffe097a20a3f5df7606375b82a4a6b..43f06735024ee732d6dd45bd24c70e40317bd74a 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1629,6 +1629,17 @@ void engine_prepare(struct engine *e) {
engine_fof(e, /*dump_results=*/0, /*seed_black_holes=*/1);
}
+ /* Perform particle splitting. Only if there is a rebuild coming and no
+ repartitioning. */
+ if (e->forcerebuild && !e->forcerepart && e->step > 1) {
+
+ /* Let's start by drifting everybody to the current time */
+ if (!drifted_all) engine_drift_all(e, /*drift_mpole=*/0);
+ drifted_all = 1;
+
+ engine_split_gas_particles(e);
+ }
+
/* Do we need repartitioning ? */
if (e->forcerepart) {
diff --git a/src/engine.h b/src/engine.h
index ca10d34467f07fa86fbb1544cc61606c655c6dfe..acd839b47f3d0fa6344b0de0f493b6e67525e58c 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -561,6 +561,9 @@ void engine_make_fof_tasks(struct engine *e);
/* Function prototypes, engine_marktasks.c. */
int engine_marktasks(struct engine *e);
+/* Function prototypes, engine_split_particles.c. */
+void engine_split_gas_particles(struct engine *e);
+
#ifdef HAVE_SETAFFINITY
cpu_set_t *engine_entry_affinity(void);
#endif
diff --git a/src/engine_split_particles.c b/src/engine_split_particles.c
new file mode 100644
index 0000000000000000000000000000000000000000..65ceb1d42d18f053b1f71206f042e80777e9de1c
--- /dev/null
+++ b/src/engine_split_particles.c
@@ -0,0 +1,296 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "engine.h"
+
+/* Local headers */
+#include "active.h"
+#include "black_holes_struct.h"
+#include "chemistry.h"
+#include "cooling.h"
+#include "hydro.h"
+#include "random.h"
+#include "star_formation.h"
+#include "tracers.h"
+
+const int particle_split_factor = 2;
+const double displacement_factor = 0.2;
+
+/**
+ * @brief Identify all the gas particles above a given mass threshold
+ * and split them into 2.
+ *
+ * This may reallocate the space arrays as new particles are created.
+ * This is an expensive function as it has to loop multiple times over
+ * the local array of particles. In case of reallocations, it may
+ * also have to loop over the gravity and other arrays.
+ *
+ * @param e The #engine.
+ */
+void engine_split_gas_particles(struct engine *e) {
+
+ /* Abort if we are not doing any splitting */
+ if (!e->hydro_properties->particle_splitting) return;
+
+ /* Time this */
+ const ticks tic = getticks();
+
+ /* Get useful constants */
+ struct space *s = e->s;
+ const int with_gravity = (e->policy & engine_policy_self_gravity) ||
+ (e->policy & engine_policy_external_gravity);
+ const double mass_threshold =
+ e->hydro_properties->particle_splitting_mass_threshold;
+
+ /* Quick check to avoid problems */
+ if (particle_split_factor != 2) {
+ error(
+ "Invalid splitting factor. Can currently only split particles into 2!");
+ }
+
+ /* Start by counting how many particles are above the threshold for splitting
+ */
+
+ size_t counter = 0;
+ const size_t nr_parts_old = s->nr_parts;
+ const struct part *parts_old = s->parts;
+ for (size_t i = 0; i < nr_parts_old; ++i) {
+
+ const struct part *p = &parts_old[i];
+
+ /* Ignore inhibited particles */
+ if (part_is_inhibited(p, e)) continue;
+
+ /* Is the mass of this particle larger than the threshold? */
+ const float gas_mass = hydro_get_mass(p);
+ if (gas_mass > mass_threshold) ++counter;
+ }
+
+ /* Early abort? */
+ if (counter == 0) {
+
+ if (e->verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+ return;
+ }
+
+ /* Be verbose about this. This is an important event */
+ message("Splitting %zd particles above the mass threshold", counter);
+
+ /* Number of particles to create */
+ const size_t count_new_gas = counter * particle_split_factor;
+
+ /* Do we need to reallocate the gas array for the new particles? */
+ if (s->nr_parts + count_new_gas > s->size_parts) {
+
+ const size_t nr_parts_new = s->nr_parts + count_new_gas;
+ s->size_parts = engine_parts_size_grow * nr_parts_new;
+
+ if (e->verbose) message("Reallocating the part array!");
+
+ struct part *parts_new = NULL;
+ if (swift_memalign("parts", (void **)&parts_new, part_align,
+ sizeof(struct part) * s->size_parts) != 0)
+ error("Failed to allocate new part data.");
+ memcpy(parts_new, s->parts, sizeof(struct part) * s->nr_parts);
+ swift_free("parts", s->parts);
+
+ struct xpart *xparts_new = NULL;
+ if (swift_memalign("xparts", (void **)&xparts_new, xpart_align,
+ sizeof(struct xpart) * s->size_parts) != 0)
+ error("Failed to allocate new part data.");
+ memcpy(xparts_new, s->xparts, sizeof(struct xpart) * s->nr_parts);
+ swift_free("xparts", s->xparts);
+
+ s->xparts = xparts_new;
+ s->parts = parts_new;
+ }
+
+ /* Do we need to reallocate the gpart array for the new particles? */
+ if (with_gravity && s->nr_gparts + count_new_gas > s->size_gparts) {
+
+ const size_t nr_gparts_new = s->nr_gparts + count_new_gas;
+ s->size_gparts = engine_parts_size_grow * nr_gparts_new;
+
+ if (e->verbose) message("Reallocating the gpart array!");
+
+ struct gpart *gparts_new = NULL;
+ if (swift_memalign("gparts", (void **)&gparts_new, gpart_align,
+ sizeof(struct gpart) * s->size_gparts) != 0)
+ error("Failed to allocate new gpart data.");
+
+ /* Offset of the new array */
+ const ptrdiff_t offset = gparts_new - s->gparts;
+
+ /* Copy the particles */
+ memcpy(gparts_new, s->gparts, sizeof(struct gpart) * s->nr_gparts);
+ swift_free("gparts", s->gparts);
+
+ /* We now need to correct all the pointers */
+ for (size_t i = 0; i < s->nr_parts; ++i) {
+ if (s->parts[i].time_bin <= num_time_bins) s->parts[i].gpart += offset;
+ }
+ /* We now need to correct all the pointers */
+ for (size_t i = 0; i < s->nr_sparts; ++i) {
+ if (s->sparts[i].time_bin <= num_time_bins) s->sparts[i].gpart += offset;
+ }
+ /* We now need to correct all the pointers */
+ for (size_t i = 0; i < s->nr_bparts; ++i) {
+ if (s->bparts[i].time_bin <= num_time_bins) s->bparts[i].gpart += offset;
+ }
+
+ s->gparts = gparts_new;
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Verify that whatever reallocation happened we are still having correct
+ * links */
+ part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, s->nr_parts,
+ s->nr_gparts, s->nr_sparts, s->nr_bparts, e->verbose);
+#endif
+
+ size_t k_parts = s->nr_parts;
+ size_t k_gparts = s->nr_gparts;
+
+ /* Loop over the particles again to split them */
+ struct part *parts = s->parts;
+ struct xpart *xparts = s->xparts;
+ struct gpart *gparts = s->gparts;
+ for (size_t i = 0; i < nr_parts_old; ++i) {
+
+ /* Is the mass of this particle larger than the threshold? */
+ struct part *p = &parts[i];
+
+ /* Ignore inhibited particles */
+ if (part_is_inhibited(p, e)) continue;
+
+ const float gas_mass = hydro_get_mass(p);
+ const float h = p->h;
+
+ /* Found a particle to split */
+ if (gas_mass > mass_threshold) {
+
+ struct xpart *xp = &xparts[i];
+ struct gpart *gp = p->gpart;
+
+ /* Start by copying over the particles */
+ memcpy(&parts[k_parts], p, sizeof(struct part));
+ memcpy(&xparts[k_parts], xp, sizeof(struct xpart));
+
+ if (with_gravity) {
+ memcpy(&gparts[k_gparts], gp, sizeof(struct gpart));
+ }
+
+ /* Update the IDs */
+ parts[k_parts].id += 2;
+
+ /* Re-link everything */
+ if (with_gravity) {
+ parts[k_parts].gpart = &gparts[k_gparts];
+ gparts[k_gparts].id_or_neg_offset = -k_parts;
+ }
+
+ /* Displacement unit vector */
+ const double delta_x = random_unit_interval(p->id, e->ti_current,
+ (enum random_number_type)0);
+ const double delta_y = random_unit_interval(p->id, e->ti_current,
+ (enum random_number_type)1);
+ const double delta_z = random_unit_interval(p->id, e->ti_current,
+ (enum random_number_type)2);
+
+ /* Displace the old particle */
+ p->x[0] += delta_x * displacement_factor * h;
+ p->x[1] += delta_y * displacement_factor * h;
+ p->x[2] += delta_z * displacement_factor * h;
+
+ if (with_gravity) {
+ gp->x[0] += delta_x * displacement_factor * h;
+ gp->x[1] += delta_y * displacement_factor * h;
+ gp->x[2] += delta_z * displacement_factor * h;
+ }
+
+ /* Displace the new particle */
+ parts[k_parts].x[0] -= delta_x * displacement_factor * h;
+ parts[k_parts].x[1] -= delta_y * displacement_factor * h;
+ parts[k_parts].x[2] -= delta_z * displacement_factor * h;
+
+ if (with_gravity) {
+ gparts[k_gparts].x[0] -= delta_x * displacement_factor * h;
+ gparts[k_gparts].x[1] -= delta_y * displacement_factor * h;
+ gparts[k_gparts].x[2] -= delta_z * displacement_factor * h;
+ }
+
+ /* Divide the mass */
+ const double new_mass = gas_mass * 0.5;
+ hydro_set_mass(p, new_mass);
+ hydro_set_mass(&parts[k_parts], new_mass);
+
+ if (with_gravity) {
+ gp->mass = new_mass;
+ gparts[k_gparts].mass = new_mass;
+ }
+
+ /* Split the chemistry fields */
+ chemistry_split_part(p, particle_split_factor);
+ chemistry_split_part(&parts[k_parts], particle_split_factor);
+
+ /* Split the cooling fields */
+ cooling_split_part(p, xp, particle_split_factor);
+ cooling_split_part(&parts[k_parts], &xparts[k_parts],
+ particle_split_factor);
+
+ /* Split the star formation fields */
+ star_formation_split_part(p, xp, particle_split_factor);
+ star_formation_split_part(&parts[k_parts], &xparts[k_parts],
+ particle_split_factor);
+
+ /* Split the tracer fields */
+ tracers_split_part(p, xp, particle_split_factor);
+ tracers_split_part(&parts[k_parts], &xparts[k_parts],
+ particle_split_factor);
+
+ /* Mark the particles as not having been swallowed */
+ black_holes_mark_part_as_not_swallowed(&p->black_holes_data);
+ black_holes_mark_part_as_not_swallowed(&parts[k_parts].black_holes_data);
+
+ /* Move to the next free slot */
+ k_parts++;
+ if (with_gravity) k_gparts++;
+ }
+ }
+
+ /* Update the local counters */
+ s->nr_parts = k_parts;
+ s->nr_gparts = k_gparts;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (s->nr_parts != nr_parts_old + (particle_split_factor - 1) * counter) {
+ error("Incorrect number of particles created");
+ }
+#endif
+
+ if (e->verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 182b8c4a803a1baf39173d08aa6d3ddc792f85da..12549d47fdfb909e4fcbfd50f6abac6579e17942 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -59,6 +59,8 @@ void hydro_props_init(struct hydro_props *p,
const struct unit_system *us,
struct swift_params *params) {
+ /* ------ Smoothing lengths parameters ---------- */
+
/* Kernel properties */
p->eta_neighbours = parser_get_param_float(params, "SPH:resolution_eta");
@@ -98,6 +100,8 @@ void hydro_props_init(struct hydro_props *p,
if (p->max_smoothing_iterations <= 10)
error("The number of smoothing length iterations should be > 10");
+ /* ------ Neighbour number definition ------------ */
+
/* Non-conventional neighbour number definition */
p->use_mass_weighted_num_ngb =
parser_get_opt_param_int(params, "SPH:use_mass_weighted_num_ngb", 0);
@@ -108,6 +112,8 @@ void hydro_props_init(struct hydro_props *p,
#endif
}
+ /* ------ Time integration parameters ------------ */
+
/* Time integration properties */
p->CFL_condition = parser_get_param_float(params, "SPH:CFL_condition");
const float max_volume_change = parser_get_opt_param_float(
@@ -121,6 +127,8 @@ void hydro_props_init(struct hydro_props *p,
if (p->initial_temperature < 0.f)
error("ERROR: Initial temperature set to a negative value!!!");
+ /* ------ Temperature parameters ----------------- */
+
/* Minimal temperature */
p->minimal_temperature = parser_get_opt_param_float(
params, "SPH:minimal_temperature", hydro_props_default_min_temp);
@@ -185,6 +193,17 @@ void hydro_props_init(struct hydro_props *p,
mean_molecular_weight = 4. / (1. + 3. * p->hydrogen_mass_fraction);
p->minimal_internal_energy = u_min / mean_molecular_weight;
+
+ /* ------ Particle splitting parameters ---------- */
+
+ /* Are we doing particle splitting? */
+ p->particle_splitting =
+ parser_get_opt_param_int(params, "SPH:particle_splitting", 0);
+
+ if (p->particle_splitting) {
+ p->particle_splitting_mass_threshold =
+ parser_get_param_float(params, "SPH:particle_splitting_mass_threshold");
+ }
}
/**
@@ -238,6 +257,12 @@ void hydro_props_print(const struct hydro_props *p) {
if (p->minimal_temperature != hydro_props_default_min_temp)
message("Minimal gas temperature set to %f", p->minimal_temperature);
+ if (p->particle_splitting)
+ message("Splitting particles with mass > %e U_M",
+ p->particle_splitting_mass_threshold);
+ else
+ message("No particle splitting");
+
/* Print out the implementation-dependent viscosity parameters
* (see hydro/SCHEME/hydro_parameters.h for this implementation) */
viscosity_print(&(p->viscosity));
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
index 54949eb91c03ed80e55d0a5a7156214b6e0e01c7..daba471b884f7d9fd24d1d85a086c59ec0431196 100644
--- a/src/hydro_properties.h
+++ b/src/hydro_properties.h
@@ -47,6 +47,8 @@ struct unit_system;
*/
struct hydro_props {
+ /* ------ Smoothing lengths parameters ---------- */
+
/*! Resolution parameter */
float eta_neighbours;
@@ -59,37 +61,43 @@ struct hydro_props {
/*! Tolerance on neighbour number (for info only)*/
float delta_neighbours;
- /*! Maximal smoothing length */
+ /*! Maximal smoothing length (internal units) */
float h_max;
/*! Minimal smoothing length expressed as ratio to softening length */
float h_min_ratio;
- /*! Minimal smoothing length */
+ /*! Minimal smoothing length (internal units) */
float h_min;
/*! Maximal number of iterations to converge h */
int max_smoothing_iterations;
+ /* ------ Neighbour number definition ------------ */
+
/*! Are we using the mass-weighted definition of neighbour number? */
int use_mass_weighted_num_ngb;
+ /* ------ Time integration parameters ------------ */
+
/*! Time integration properties */
float CFL_condition;
/*! Maximal change of h over one time-step */
float log_max_h_change;
+ /* ------ Temperature parameters ----------------- */
+
/*! Minimal temperature allowed */
float minimal_temperature;
- /*! Minimal physical internal energy per unit mass */
+ /*! Minimal physical internal energy per unit mass (internal units) */
float minimal_internal_energy;
/*! Initial temperature */
float initial_temperature;
- /*! Initial physical internal energy per unit mass */
+ /*! Initial physical internal energy per unit mass (internal units) */
float initial_internal_energy;
/*! Primordial hydrogen mass fraction for initial energy conversion */
@@ -104,6 +112,16 @@ struct hydro_props {
/*! Mean molecular weight above hydrogen ionization temperature */
float mu_ionised;
+ /* ------ Particle splitting parameters ---------- */
+
+ /*! Is particle splitting activated? */
+ int particle_splitting;
+
+ /*! Mass above which particles get split (internal units) */
+ float particle_splitting_mass_threshold;
+
+ /* ------ Viscosity and diffusion ---------------- */
+
/*! Artificial viscosity parameters */
struct viscosity_global_data viscosity;
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index 767fa2715ab35f21ae41c5162787c3a544b87da7..143984ebf8758c2510b4f1eca9f0a1e6b912896e 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/src/star_formation/EAGLE/star_formation.h
@@ -722,4 +722,20 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const,
const struct part* restrict p,
struct xpart* restrict xp) {}
+/**
+ * @brief Split the star formation content of a particle into n pieces
+ *
+ * We only need to split the SFR if it is positive, i.e. it is not
+ * storing the redshift/time of last SF event.
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+__attribute__((always_inline)) INLINE static void star_formation_split_part(
+ struct part* p, struct xpart* xp, const double n) {
+
+ if (xp->sf_data.SFR > 0.) xp->sf_data.SFR /= n;
+}
+
#endif /* SWIFT_EAGLE_STAR_FORMATION_H */
diff --git a/src/star_formation/EAGLE/star_formation_struct.h b/src/star_formation/EAGLE/star_formation_struct.h
index 41247e160a3eddbc9184c59b67cfa2a1d7259a05..f5f0f76e8d59f73148f51f462cc619ec70e82192 100644
--- a/src/star_formation/EAGLE/star_formation_struct.h
+++ b/src/star_formation/EAGLE/star_formation_struct.h
@@ -25,7 +25,8 @@
*/
struct star_formation_xpart_data {
- /*! Star formation rate */
+ /*! Star formation rate (internal units) or (if negative) time/scale-factor of
+ * last SF episode */
float SFR;
};
diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h
index 1f343f5e649c557ad6bdff7e9fcb377eea796f29..c2682f1643f6a998e5e9ec1a58b5f85d9484af51 100644
--- a/src/star_formation/GEAR/star_formation.h
+++ b/src/star_formation/GEAR/star_formation.h
@@ -287,4 +287,17 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const,
const struct part* restrict p,
struct xpart* restrict xp) {}
+/**
+ * @brief Split the star formation content of a particle into n pieces
+ *
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+__attribute__((always_inline)) INLINE static void star_formation_split_part(
+ struct part* p, struct xpart* xp, const double n) {
+ error("Loic: to be implemented");
+}
+
#endif /* SWIFT_GEAR_STAR_FORMATION_H */
diff --git a/src/star_formation/none/star_formation.h b/src/star_formation/none/star_formation.h
index 640a6c378839bed8c5ac824d0cc77773cb4e728f..a48155b05beb7a345a006ebeb09961a66483ba40 100644
--- a/src/star_formation/none/star_formation.h
+++ b/src/star_formation/none/star_formation.h
@@ -227,4 +227,14 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const,
const struct part* restrict p,
struct xpart* restrict xp) {}
+/**
+ * @brief Split the star formation content of a particle into n pieces
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+__attribute__((always_inline)) INLINE static void star_formation_split_part(
+ struct part* p, struct xpart* xp, const double n) {}
+
#endif /* SWIFT_NONE_STAR_FORMATION_H */
diff --git a/src/tracers/EAGLE/tracers.h b/src/tracers/EAGLE/tracers.h
index 046a2088fef3f595e5ae4ab72b66f54a2dca991c..57947b37f0c3c85a37ca0d6e24a9575d58348b80 100644
--- a/src/tracers/EAGLE/tracers.h
+++ b/src/tracers/EAGLE/tracers.h
@@ -158,4 +158,16 @@ static INLINE void tracers_after_black_holes_feedback(struct xpart *xp) {
xp->tracers_data.hit_by_AGN_feedback = 1;
}
+/**
+ * @brief Split the tracer content of a particle into n pieces
+ *
+ * Nothing to do here.
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+__attribute__((always_inline)) INLINE static void tracers_split_part(
+ struct part *p, struct xpart *xp, const double n) {}
+
#endif /* SWIFT_TRACERS_EAGLE_H */
diff --git a/src/tracers/none/tracers.h b/src/tracers/none/tracers.h
index affb9c255a126333c7bd5cfc5ed2b42318652169..99bf395b92dbc5185c2aee1e26f41c38971797c7 100644
--- a/src/tracers/none/tracers.h
+++ b/src/tracers/none/tracers.h
@@ -127,4 +127,15 @@ static INLINE void tracers_after_feedback(struct xpart *xp) {}
*/
static INLINE void tracers_after_black_holes_feedback(struct xpart *xp) {}
+/**
+ * @brief Split the tracer content of a particle into n pieces
+ *
+ * Nothing to do here.
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param n The number of pieces to split into.
+ */
+__attribute__((always_inline)) INLINE static void tracers_split_part(
+ struct part *p, struct xpart *xp, const double n) {}
#endif /* SWIFT_TRACERS_NONE_H */
diff --git a/tools/analyse_runtime.py b/tools/analyse_runtime.py
index e94d10e17ad652dd49a82d3562b9a8ea31670df0..f0fee27266d0d20ae9234090581676a9c8df6860 100755
--- a/tools/analyse_runtime.py
+++ b/tools/analyse_runtime.py
@@ -54,6 +54,7 @@ threshold = 0.008
num_files = len(sys.argv) - 1
labels = [
+ ["engine_split_gas_particles:", 1],
["Gpart assignment", 1],
["Mesh comunication", 1],
["Forward Fourier transform", 1],