diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/makeDisk.py b/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/makeDisk.py
index cec23105e59edc4e1a6c551a80b722b4eb580d05..70f5c2d9c49293210f185b936b4d264ba82fd462 100755
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/makeDisk.py
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/makeDisk.py
@@ -588,14 +588,14 @@ nb.rename("galaxy_multi_component.hdf5")
nb.write()
#%%
-#Add the StellarParticleType attribute to the dataset
+# Add the StellarParticleType attribute to the dataset
import h5py as h5
import numpy as np
nb_star = nb.select("stars")
N_star = np.sum(nb_star.npart)
-star_tpe = 2 # Single population stars
-star_type = np.ones(N_star)*star_tpe
+star_tpe = 2 # Single population stars
+star_type = np.ones(N_star) * star_tpe
-with h5.File("galaxy_multi_component.hdf5", "r+") as f:
+with h5.File("galaxy_multi_component.hdf5", "r+") as f:
f["PartType4"].create_dataset("StellarParticleType", data=star_type)
diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/README b/examples/SinkParticles/HomogeneousBoxSinkParticles/README
new file mode 100644
index 0000000000000000000000000000000000000000..8669089494c529ed3e2fbe4c9ebc8c084574eeea
--- /dev/null
+++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/README
@@ -0,0 +1,24 @@
+# Intro
+This example is a modified version of the HomogeneousBox in SWIFT that allows to add sink particles to the box. By default, there is no gravity.
+
+The default level is 5 (N_particle = 2^(3*level)). It is quick. If you want to try the example with higher levels, we recommend using HPC facilities.
+
+*This example is mainly used for sink MPI debugging.*
+
+# Configure
+
+To run this example with GEAR model,
+
+./configure --with-chemistry=GEAR_10 --with-cooling=grackle_0 --with-stars=GEAR --with-star-formation=GEAR --with-feedback=GEAR --with-sink=GEAR --with-kernel=wendland-C2 --with-grackle=path/to/grackle --enable-debug --enable-debugging-checks
+
+and then
+
+make -j
+
+# ICs
+The run.sh script calls `makeIC.py' script with default values. You can experiment by changing the ICs. Run `python3 makeIC.py --help` to get the list of parameters.
+
+# Run
+Type `run.sh` (or `n_ranks=4 mpi_run.sh`), and let's go!
+
+You can provide parameters to the running scripts, have a look inside them.
diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/getChemistryTable.sh b/examples/SinkParticles/HomogeneousBoxSinkParticles/getChemistryTable.sh
new file mode 100755
index 0000000000000000000000000000000000000000..b10fd23964158ee7a38d352dbd0ddd9beb7bdd77
--- /dev/null
+++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/getChemistryTable.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/FeedbackTables/POPIIsw.h5
diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/getGrackleCoolingTable.sh b/examples/SinkParticles/HomogeneousBoxSinkParticles/getGrackleCoolingTable.sh
new file mode 100755
index 0000000000000000000000000000000000000000..e3eb106240709c80151a48625567d2cd78e5f568
--- /dev/null
+++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/getGrackleCoolingTable.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5
diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py b/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..6ad3be507b4fd12d8bef3c0af4a1b674cb950a27
--- /dev/null
+++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py
@@ -0,0 +1,267 @@
+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Yves Revaz (yves.revaz@epfl.ch)
+# 2024 Darwin Roduit (darwin.roduit@epfl.ch)
+#
+# 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/>.
+#
+################################################################################
+
+import h5py
+import numpy as np
+import argparse
+from astropy import units
+from astropy import constants
+
+
+class store_as_array(argparse._StoreAction):
+ """Provides numpy array as argparse arguments."""
+
+ def __call__(self, parser, namespace, values, option_string=None):
+ values = np.array(values)
+ return super().__call__(parser, namespace, values, option_string)
+
+
+def parse_options():
+
+ usage = "usage: %prog [options] file"
+ parser = argparse.ArgumentParser(description=usage)
+
+ parser.add_argument(
+ "--lJ",
+ action="store",
+ dest="lJ",
+ type=float,
+ default=0.250,
+ help="Jeans wavelength in box size unit",
+ )
+
+ parser.add_argument(
+ "--rho",
+ action="store",
+ dest="rho",
+ type=float,
+ default=0.1,
+ help="Mean gas density in atom/cm3",
+ )
+
+ parser.add_argument(
+ "--mass",
+ action="store",
+ dest="mass",
+ type=float,
+ default=50,
+ help="Gas particle mass in solar mass",
+ )
+
+ parser.add_argument(
+ "--level",
+ action="store",
+ dest="level",
+ type=int,
+ default=5,
+ help="Resolution level: N = (2**l)**3",
+ )
+
+ parser.add_argument(
+ "--sink_mass",
+ action="store",
+ type=float,
+ default=50,
+ help="Sink particles mass in solar mass",
+ )
+
+ parser.add_argument(
+ "--sinks_vel",
+ action=store_as_array,
+ nargs=3,
+ type=float,
+ default=np.array([10, 0, 0]),
+ help="Sink particle velocity. All sinks get the same velocity",
+ )
+
+ parser.add_argument(
+ "--sink_pos",
+ action=store_as_array,
+ nargs=3,
+ type=float,
+ default=np.array([0, 0, 0]),
+ help="Sink particle position. Only use it to place one sink particle.",
+ )
+
+ parser.add_argument(
+ "--n_sink",
+ action="store",
+ type=int,
+ default=10,
+ help="Number of sink particles in the box",
+ )
+
+ parser.add_argument(
+ "-o",
+ action="store",
+ dest="outputfilename",
+ type=str,
+ default="box.hdf5",
+ help="output filename",
+ )
+
+ options = parser.parse_args()
+ return options
+
+
+########################################
+# main
+########################################
+
+opt = parse_options()
+
+# define standard units
+UnitMass_in_cgs = 1.988409870698051e43 # 10^10 M_sun in grams
+UnitLength_in_cgs = 3.0856775814913673e21 # kpc in centimeters
+UnitVelocity_in_cgs = 1e5 # km/s in centimeters per second
+UnitCurrent_in_cgs = 1 # Amperes
+UnitTemp_in_cgs = 1 # Kelvin
+UnitTime_in_cgs = UnitLength_in_cgs / UnitVelocity_in_cgs
+
+UnitMass = UnitMass_in_cgs * units.g
+UnitLength = UnitLength_in_cgs * units.cm
+UnitTime = UnitTime_in_cgs * units.s
+UnitVelocity = UnitVelocity_in_cgs * units.cm / units.s
+
+np.random.seed(1)
+
+# Number of particles
+N = (2 ** opt.level) ** 3 # number of particles
+
+# Mean density
+rho = opt.rho # atom/cc
+rho = rho * constants.m_p / units.cm ** 3
+
+# Gas particle mass
+m = opt.mass # in solar mass
+m = m * units.Msun
+
+# Gas mass in the box
+M = N * m
+
+# Size of the box
+L = (M / rho) ** (1 / 3.0)
+
+# Jeans wavelength in box size unit
+lJ = opt.lJ
+lJ = lJ * L
+
+# Gravitational constant
+G = constants.G
+
+# Jeans wave number
+kJ = 2 * np.pi / lJ
+# Velocity dispersion
+sigma = np.sqrt(4 * np.pi * G * rho) / kJ
+
+
+print("Boxsize : {}".format(L.to(units.kpc)))
+print("Number of particles : {}".format(N))
+print("Equivalent velocity dispertion : {}".format(sigma.to(units.m / units.s)))
+
+# Convert to code units
+m = m.to(UnitMass).value
+L = L.to(UnitLength).value
+rho = rho.to(UnitMass / UnitLength ** 3).value
+sigma = sigma.to(UnitVelocity).value
+
+# Generate the particles
+pos = np.random.random([N, 3]) * np.array([L, L, L])
+vel = np.zeros([N, 3])
+mass = np.ones(N) * m
+u = np.ones(N) * sigma ** 2
+ids = np.arange(N)
+h = np.ones(N) * 3 * L / N ** (1 / 3.0)
+rho = np.ones(N) * rho
+
+print("Inter-particle distance (code unit) : {}".format(L / N ** (1 / 3.0)))
+
+
+#####################
+# Now, take care of the sink
+#####################
+N_sink = opt.n_sink
+m_sink = opt.sink_mass * units.M_sun
+m_sink = m_sink.to(UnitMass).value # Convert the sink mass to internal units
+
+
+if N_sink == 1:
+ pos_sink = np.reshape(opt.sinks_vel, (N_sink, 3))
+else:
+ pos_sink = np.random.random([N_sink, 3]) * np.array([L, L, L])
+
+if opt.sinks_vel is not None:
+ vel_sink = np.tile(
+ opt.sinks_vel, (N_sink, 1)
+ ) # Duplicate the velocity for all sinks
+else:
+ np.zeros([N_sink, 3])
+
+mass_sink = np.ones(N_sink) * m_sink
+ids_sink = np.arange(N, N + N_sink)
+
+#####################
+# Finally write the ICs in the file
+#####################
+
+# File
+fileOutput = h5py.File(opt.outputfilename, "w")
+print("{} saved.".format(opt.outputfilename))
+
+# Header
+grp = fileOutput.create_group("/Header")
+grp.attrs["BoxSize"] = [L, L, L]
+grp.attrs["NumPart_Total"] = [N, 0, 0, N_sink, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [N, 0, 0, N_sink, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFileOutputsPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["Dimension"] = 3
+
+
+# Units
+grp = fileOutput.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = UnitLength_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = UnitMass_in_cgs
+grp.attrs["Unit time in cgs (U_t)"] = UnitTime_in_cgs
+grp.attrs["Unit current in cgs (U_I)"] = UnitCurrent_in_cgs
+grp.attrs["Unit temperature in cgs (U_T)"] = UnitTemp_in_cgs
+
+
+# Write Gas particle group
+grp = fileOutput.create_group("/PartType0")
+grp.create_dataset("Coordinates", data=pos, dtype="d")
+grp.create_dataset("Velocities", data=vel, dtype="f")
+grp.create_dataset("Masses", data=mass, dtype="f")
+grp.create_dataset("SmoothingLength", data=h, dtype="f")
+grp.create_dataset("InternalEnergy", data=u, dtype="f")
+grp.create_dataset("ParticleIDs", data=ids, dtype="L")
+grp.create_dataset("Densities", data=rho, dtype="f")
+
+
+# Write sink particle group
+grp = fileOutput.create_group("/PartType3")
+grp.create_dataset("Coordinates", data=pos_sink, dtype="d")
+grp.create_dataset("Velocities", data=vel_sink, dtype="f")
+grp.create_dataset("Masses", data=mass_sink, dtype="f")
+grp.create_dataset("ParticleIDs", data=ids_sink, dtype="L")
+fileOutput.close()
diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/mpi_run.sh b/examples/SinkParticles/HomogeneousBoxSinkParticles/mpi_run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..8b884691d68c92902a68b62826c4a7dcfaa9dc91
--- /dev/null
+++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/mpi_run.sh
@@ -0,0 +1,51 @@
+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+
+n_ranks=${n_ranks:=2} #Number of MPI ranks
+n_threads=${n_threads:=1} #Number of threads to use
+level=${level:=5} #Number of particles = 2^(3*level)
+jeans_length=${jeans_length:=0.250} #Jeans wavelenght in unit of the boxsize
+n_sinks=${n_sinks:=10} #Number of sinks
+
+# Remove the ICs
+if [ -e ICs_homogeneous_box.hdf5 ]
+then
+ rm ICs_homogeneous_box.hdf5
+fi
+
+#Create the ICs if they do not exist
+if [ ! -e ICs_homogeneous_box.hdf5 ]
+then
+ echo "Generating initial conditions to run the example..."
+ python3 makeIC.py --level $level -o ICs_homogeneous_box.hdf5 --lJ $jeans_length --n_sink $n_sinks --sink_pos 0 0 0 --sinks_vel 10 10 0
+fi
+
+# Get the Grackle cooling table
+if [ ! -e CloudyData_UVB=HM2012.h5 ]
+then
+ echo "Fetching the Cloudy tables required by Grackle..."
+ ./getGrackleCoolingTable.sh
+fi
+
+if [ ! -e POPIIsw.h5 ]
+then
+ echo "Fetching the chemistry tables..."
+ ./getChemistryTable.sh
+fi
+
+# Create output directory
+DIR=snap #First test of units conversion
+if [ -d "$DIR" ];
+then
+ echo "$DIR directory exists. Its content will be removed."
+ rm -r $DIR
+else
+ echo "$DIR directory does not exists. It will be created."
+ mkdir $DIR
+fi
+
+printf "Running simulation..."
+
+mpirun -n $n_ranks ../../../swift_mpi --pin --hydro --sinks --stars --external-gravity --feedback --threads=$n_threads --verbose 1 params.yml 2>&1 | tee output.log
diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/params.yml b/examples/SinkParticles/HomogeneousBoxSinkParticles/params.yml
new file mode 100755
index 0000000000000000000000000000000000000000..fd7c945e8dadf6ba204f973f60d652981cfa4bee
--- /dev/null
+++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/params.yml
@@ -0,0 +1,101 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.988409870698051e+43 # 10^10 solar masses
+ UnitLength_in_cgs: 3.0856775814913673e+21 # 1 kpc
+ UnitVelocity_in_cgs: 1e5 # km/s
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters for the self-gravity scheme
+Gravity:
+ MAC: adaptive # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'.
+ epsilon_fmm: 0.001 # Tolerance parameter for the adaptive multipole acceptance criterion.
+ theta_cr: 0.7 # Opening angle for the purely gemoetric criterion.
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.7 # Opening angle (Multipole acceptance criterion).
+ max_physical_baryon_softening: 0.005 # Physical softening length (in internal units)
+ mesh_side_length: 32
+
+# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 0.282 # 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: 2e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ subdir: snap
+ basename: snapshot # Common part of the name of output files
+ time_first: 0. #230 Myr # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
+ delta_time: 10e-3 # Time difference between consecutive outputs (in internal units)
+
+Scheduler:
+ cell_extra_gparts: 100 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
+ cell_extra_sinks: 100 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
+ cell_extra_sparts: 100 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
+ max_top_level_cells: 3 #
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 5e-3 # Time between statistics output
+ time_first: 0. # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units)
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ICs_homogeneous_box.hdf5
+ periodic: 1 # Are we running with periodic ICs?
+ shift: [0.0,0.0,0.0]
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 57Ngbs with the Wendland C2 kernel).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ h_max: 5
+ minimal_temperature: 1
+
+
+# Cooling with Grackle 3.0
+GrackleCooling:
+ cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+ with_UV_background: 0 # Enable or not the UV background
+ redshift: -1 # Redshift to use (-1 means time based redshift)
+ with_metal_cooling: 1 # Enable or not the metal cooling
+ provide_volumetric_heating_rates: 0 # User provide volumetric heating rates
+ provide_specific_heating_rates: 0 # User provide specific heating rates
+ self_shielding_method: -1 # Grackle (<= 3) or Gear self shielding method
+ self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
+ max_steps: 1000
+ convergence_limit: 1e-2
+ thermal_time_myr: 5
+ maximal_density_Hpcm3: -1 # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter.
+
+GEARChemistry:
+ initial_metallicity: 0
+
+GEARFeedback:
+ supernovae_energy_erg: 1e51 # supernovae energy, used only for SNIa
+ supernovae_efficiency: 0.1 # supernovae energy efficiency, used for both SNIa and SNII
+ yields_table: POPIIsw.h5
+ yields_table_first_stars: POPIIsw.h5
+ discrete_yields: 1
+ imf_transition_metallicity: -5 # Maximal metallicity ([Fe/H]) for a first star (0 to deactivate).
+ elements: [Fe, Mg, O, C, Al, Ca, Ba, Zn, Eu] # Elements to read in the yields table. The number of element should be one less than the number of elements (N) requested during the configuration (--with-chemistry=GEAR_N).
+
+GEARSink:
+ cut_off_radius: 5e-3 # Cut off radius of all the sinks in internal units.
+ f_acc: 0.1
+ temperature_threshold_K: 100 # Max temperature (in K) for forming a sink when density_threshold_g_per_cm3 <= density <= maximal_density_threshold_g_per_cm3.
+ density_threshold_Hpcm3: 1e1 # Minimum gas density (in g/cm3) required to form a sink particle.
+ maximal_density_threshold_Hpcm3: 1e5 # If the gas density exceeds this value (in g/cm3), a sink forms regardless of temperature if all other criteria are passed.
+ stellar_particle_mass_Msun: 50 # Mass of the stellar particle representing the low mass stars, in solar mass
+ minimal_discrete_mass_Msun: 8 # Minimal mass of stars represented by discrete particles, in solar mass
+ stellar_particle_mass_first_stars_Msun: 50 # Mass of the stellar particle representing the low mass stars, in solar mass
+ minimal_discrete_mass_first_stars_Msun: 8 # Minimal mass of stars represented by discrete particles, in solar mass
+ star_spawning_sigma_factor: 0.5 # Factor to rescale the velocity dispersion of the stars when they are spawned. (Default: 0.2)
+ sink_formation_contracting_gas_criterion: 1 # (Optional) Activate the contracting gas check for sink formation. (Default: 1)
+ sink_formation_smoothing_length_criterion: 1 # (Optional) Activate the smoothing length check for sink formation. (Default: 1)
+ sink_formation_jeans_instability_criterion: 1 # (Optional) Activate the two Jeans instability checks for sink formation. (Default: 1)
+ sink_formation_bound_state_criterion: 1 # (Optional) Activate the bound state check for sink formation. (Default: 1)
+ sink_formation_overlapping_sink_criterion: 1 # (Optional) Activate the overlapping sink check for sink formation. (Default: 1)
+ disable_sink_formation: 0 # (Optional) Disable sink formation. (Default: 0)
diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/run.sh b/examples/SinkParticles/HomogeneousBoxSinkParticles/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..9db866d75f792434d5fe2d93f9e6668299737dd2
--- /dev/null
+++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/run.sh
@@ -0,0 +1,55 @@
+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+
+n_threads=${n_threads:=8} #Number of threads to use
+level=${level:=5} #Number of particles = 2^(3*level)
+jeans_length=${jeans_length:=0.250} #Jeans wavelenght in unit of the boxsize
+gas_density=${gas_density:=0.1} #Gas density in atom/cm^3
+gas_particle_mass=${gas_particle_mass:=50} #Mass of the gas particles
+n_sinks=${n_sinks:=10} #Number of sinks
+
+# Remove the ICs
+if [ -e ICs_homogeneous_box.hdf5 ]
+then
+ rm ICs_homogeneous_box.hdf5
+fi
+
+#Create the ICs if they do not exist
+if [ ! -e ICs_homogeneous_box.hdf5 ]
+then
+ echo "Generating initial conditions to run the example..."
+ python3 makeIC.py --level $level -o ICs_homogeneous_box.hdf5 --lJ $jeans_length --n_sink $n_sinks --sink_pos 0 0 0 --sinks_vel 10 10 0
+fi
+
+
+# Get the Grackle cooling table
+if [ ! -e CloudyData_UVB=HM2012.h5 ]
+then
+ echo "Fetching the Cloudy tables required by Grackle..."
+ ./getGrackleCoolingTable.sh
+fi
+
+
+if [ ! -e POPIIsw.h5 ]
+then
+ echo "Fetching the chemistry tables..."
+ ./getChemistryTable.sh
+fi
+
+
+# Create output directory
+DIR=snap #First test of units conversion
+if [ -d "$DIR" ];
+then
+ echo "$DIR directory exists. Its content will be removed."
+ rm -r $DIR
+else
+ echo "$DIR directory does not exists. It will be created."
+ mkdir $DIR
+fi
+
+printf "Running simulation..."
+
+../../../swift --hydro --sinks --stars --external-gravity --feedback --threads=$n_threads params.yml 2>&1 | tee output.log
diff --git a/src/cell.c b/src/cell.c
index 4f58054830f31a43dd96f8539262af0fc0b5ceb6..faadb774325f5bdf38ddc1c43742ed261c84a45f 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -293,6 +293,39 @@ int cell_link_bparts(struct cell *c, struct bpart *bparts) {
return c->black_holes.count;
}
+/**
+ * @brief Link the cells recursively to the given #sink array.
+ *
+ * @param c The #cell.
+ * @param sinks The #sink array.
+ *
+ * @return The number of particles linked.
+ */
+int cell_link_sinks(struct cell *c, struct sink *sinks) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->nodeID == engine_rank)
+ error("Linking foreign particles in a local cell!");
+
+ if (c->sinks.parts != NULL)
+ error("Linking sparts into a cell that was already linked");
+#endif
+
+ c->sinks.parts = sinks;
+ c->sinks.parts_rebuild = sinks;
+
+ /* Fill the progeny recursively, depth-first. */
+ if (c->split) {
+ int offset = 0;
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL)
+ offset += cell_link_sinks(c->progeny[k], &sinks[offset]);
+ }
+ }
+
+ /* Return the total number of linked particles. */
+ return c->sinks.count;
+}
+
/**
* @brief Recurse down foreign cells until reaching one with hydro
* tasks; then trigger the linking of the #part array from that
@@ -413,6 +446,7 @@ void cell_unlink_foreign_particles(struct cell *c) {
c->hydro.parts = NULL;
c->stars.parts = NULL;
c->black_holes.parts = NULL;
+ c->sinks.parts = NULL;
if (c->split) {
for (int k = 0; k < 8; k++) {
diff --git a/src/cell.h b/src/cell.h
index 2bb2d6976964461d797069a24c9bad9f9e71cdf3..9cb04b32c4a81affa50ea6584c1679ca756414dc 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -280,6 +280,15 @@ struct pcell_step {
float dx_max_part;
} black_holes;
+ struct {
+
+ /*! Minimal integer end-of-timestep in this cell (sinks) */
+ integertime_t ti_end_min;
+
+ /*! Maximal distance any #part has travelled since last rebuild */
+ float dx_max_part;
+ } sinks;
+
struct {
/*! Minimal integer end-of-timestep in this cell (rt) */
@@ -562,6 +571,7 @@ int cell_link_parts(struct cell *c, struct part *parts);
int cell_link_gparts(struct cell *c, struct gpart *gparts);
int cell_link_sparts(struct cell *c, struct spart *sparts);
int cell_link_bparts(struct cell *c, struct bpart *bparts);
+int cell_link_sinks(struct cell *c, struct sink *sinks);
int cell_link_foreign_parts(struct cell *c, struct part *parts);
int cell_link_foreign_gparts(struct cell *c, struct gpart *gparts);
void cell_unlink_foreign_particles(struct cell *c);
diff --git a/src/cell_pack.c b/src/cell_pack.c
index 21bf9ad122c788f228e811a3dd7024d39fa5cafb..83f8c4cab9d9c9de3cddfc906f789c3e0d2b2384 100644
--- a/src/cell_pack.c
+++ b/src/cell_pack.c
@@ -303,6 +303,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
temp->hydro.count = 0;
temp->grav.count = 0;
temp->stars.count = 0;
+ temp->sinks.count = 0;
temp->loc[0] = c->loc[0];
temp->loc[1] = c->loc[1];
temp->loc[2] = c->loc[2];
@@ -319,6 +320,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
temp->hydro.dx_max_sort = 0.f;
temp->stars.dx_max_part = 0.f;
temp->stars.dx_max_sort = 0.f;
+ temp->sinks.dx_max_part = 0.f;
temp->black_holes.dx_max_part = 0.f;
temp->nodeID = c->nodeID;
temp->parent = c;
@@ -457,6 +459,9 @@ int cell_pack_end_step(const struct cell *c, struct pcell_step *pcells) {
pcells[0].black_holes.ti_end_min = c->black_holes.ti_end_min;
pcells[0].black_holes.dx_max_part = c->black_holes.dx_max_part;
+ pcells[0].sinks.ti_end_min = c->sinks.ti_end_min;
+ pcells[0].sinks.dx_max_part = c->sinks.dx_max_part;
+
/* Fill in the progeny, depth-first recursion. */
int count = 1;
for (int k = 0; k < 8; k++)
@@ -501,6 +506,9 @@ int cell_unpack_end_step(struct cell *c, const struct pcell_step *pcells) {
c->black_holes.ti_end_min = pcells[0].black_holes.ti_end_min;
c->black_holes.dx_max_part = pcells[0].black_holes.dx_max_part;
+ c->sinks.ti_end_min = pcells[0].sinks.ti_end_min;
+ c->sinks.dx_max_part = pcells[0].sinks.dx_max_part;
+
/* Fill in the progeny, depth-first recursion. */
int count = 1;
for (int k = 0; k < 8; k++)
diff --git a/src/cell_sinks.h b/src/cell_sinks.h
index 6bd2f9a9afed9847001a9a106b6f2dac9a1b433d..9043ad191552657c7183b35ec7d443207f4af384 100644
--- a/src/cell_sinks.h
+++ b/src/cell_sinks.h
@@ -42,6 +42,9 @@ struct cell_sinks {
/*! Pointer to the #sink data. */
struct sink *parts;
+ /*! Pointer to the #spart data at rebuild time. */
+ struct sink *parts_rebuild;
+
/*! Linked list of the tasks computing this cell's sink swallow. */
struct link *swallow;
diff --git a/src/cell_split.c b/src/cell_split.c
index ee669ae17ac05eb8de1c8724c4656d8946d04c33..df0fc2572e0b5fa1c33164c07b7b1c7453a3804e 100644
--- a/src/cell_split.c
+++ b/src/cell_split.c
@@ -384,6 +384,7 @@ void cell_split(struct cell *c, const ptrdiff_t parts_offset,
c->progeny[k]->sinks.count = bucket_count[k];
c->progeny[k]->sinks.count_total = c->progeny[k]->sinks.count;
c->progeny[k]->sinks.parts = &c->sinks.parts[bucket_offset[k]];
+ c->progeny[k]->sinks.parts_rebuild = c->progeny[k]->sinks.parts;
}
/* Finally, do the same song and dance for the gparts. */
diff --git a/src/cell_unskip.c b/src/cell_unskip.c
index 3ce7d791a712e5dcacbc7794af08e81513fa5fdd..adfd6eb6ef1512c845a46674681cc6c1edea0bc3 100644
--- a/src/cell_unskip.c
+++ b/src/cell_unskip.c
@@ -3001,7 +3001,7 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
if (cell_need_rebuild_for_sinks_pair(ci, cj)) rebuild = 1;
if (cell_need_rebuild_for_sinks_pair(cj, ci)) rebuild = 1;
-#ifdef WITH_MPI
+#if defined(WITH_MPI) && !defined(SWIFT_DEBUG_CHECKS)
error("TODO");
#endif
}
diff --git a/src/cooling/grackle/cooling.c b/src/cooling/grackle/cooling.c
index 205a7b4f96ca663fbf0ee540febb38d4f5479488..be7d417259654ccdf370fb644de1db012f5f683a 100644
--- a/src/cooling/grackle/cooling.c
+++ b/src/cooling/grackle/cooling.c
@@ -68,6 +68,18 @@ double cooling_get_physical_density(
const struct part* p, const struct cosmology* cosmo,
const struct cooling_function_data* cooling);
+/**
+ * @brief Record the time when cooling was switched off for a particle.
+ *
+ * @param p #part data.
+ * @param xp Pointer to the #xpart data.
+ * @param time The time when the cooling was switched off.
+ */
+INLINE void cooling_set_part_time_cooling_off(struct part* p, struct xpart* xp,
+ const double time) {
+ xp->cooling_data.time_last_event = time;
+}
+
/**
* @brief Common operations performed on the cooling function at a
* given time-step or redshift.
diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h
index 92f6c921634bb62e58626effbcdc1019b2ecc75f..e5849a421eff0763217b09c8bcf15734471a1b94 100644
--- a/src/cooling/grackle/cooling.h
+++ b/src/cooling/grackle/cooling.h
@@ -105,6 +105,16 @@ INLINE static float cooling_get_subgrid_density(const struct part* p,
return -1.f;
}
+/**
+ * @brief Record the time when cooling was switched off for a particle.
+ *
+ * @param p #part data.
+ * @param xp Pointer to the #xpart data.
+ * @param time The time when the cooling was switched off.
+ */
+void cooling_set_part_time_cooling_off(struct part* p, struct xpart* xp,
+ const double time);
+
float cooling_get_radiated_energy(const struct xpart* restrict xp);
void cooling_print_backend(const struct cooling_function_data* cooling);
diff --git a/src/cooling/none/cooling.h b/src/cooling/none/cooling.h
index dbb5e9872bf4e197f372c3f4a121c1beaaa775f9..d7955c4e4b1902796b855012e3f454f69e0eb70e 100644
--- a/src/cooling/none/cooling.h
+++ b/src/cooling/none/cooling.h
@@ -277,6 +277,19 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
return 0.f;
}
+/**
+ * @brief Record the time when cooling was switched off for a particle.
+ *
+ * In none model, this function does nothing.
+ *
+ * @param p #part data.
+ * @param xp Pointer to the #xpart data.
+ * @param time The time when the cooling was switched off.
+ */
+INLINE static void cooling_set_part_time_cooling_off(struct part* p,
+ struct xpart* xp,
+ const double time) {}
+
/**
* @brief Split the coolong content of a particle into n pieces
*
diff --git a/src/engine.c b/src/engine.c
index 3b992186dcdf2e9907accda1598cc084f4156603..1707dcf6300f1224e9fc8ccac7b769206653c64f 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -781,13 +781,14 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
e->policy & (engine_policy_hydro | engine_policy_grid_hydro);
const int with_stars = e->policy & engine_policy_stars;
const int with_black_holes = e->policy & engine_policy_black_holes;
+ const int with_sinks = e->policy & engine_policy_sinks;
struct space *s = e->s;
ticks tic = getticks();
/* Count the number of particles we need to import and re-allocate
the buffer if needed. */
size_t count_parts_in = 0, count_gparts_in = 0, count_sparts_in = 0,
- count_bparts_in = 0;
+ count_bparts_in = 0, count_sinks_in = 0;
for (int k = 0; k < nr_proxies; k++) {
for (int j = 0; j < e->proxies[k].nr_cells_in; j++) {
@@ -806,6 +807,11 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
/* For black holes, we just use the numbers in the top-level cells */
count_bparts_in += e->proxies[k].cells_in[j]->black_holes.count;
+
+ /* For sinks, we just use the numbers in the top-level cells + some
+ extra space */
+ count_sinks_in +=
+ e->proxies[k].cells_in[j]->sinks.count + space_extra_sinks;
}
}
@@ -871,28 +877,49 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
error("Failed to allocate foreign bpart data.");
}
+ /* Allocate space for the foreign particles we will receive */
+ size_t old_size_sinks_foreign = s->size_sinks_foreign;
+ if (!fof && count_sinks_in > s->size_sinks_foreign) {
+ if (s->sinks_foreign != NULL) swift_free("sinks_foreign", s->sinks_foreign);
+ s->size_sinks_foreign = engine_foreign_alloc_margin * count_sinks_in;
+ if (swift_memalign("sinks_foreign", (void **)&s->sinks_foreign, sink_align,
+ sizeof(struct sink) * s->size_sinks_foreign) != 0)
+ error("Failed to allocate foreign sink data.");
+ bzero(s->sinks_foreign, s->size_sinks_foreign * sizeof(struct sink));
+
+ /* Note: If you ever see a sink particle with id = -666, the following
+ lines is the ones that sets the ID to this value. */
+ for (size_t i = 0; i < s->size_sinks_foreign; ++i) {
+ s->sinks_foreign[i].time_bin = time_bin_not_created;
+ s->sinks_foreign[i].id = -666;
+ }
+ }
+
if (e->verbose) {
message(
- "Allocating %zd/%zd/%zd/%zd foreign part/gpart/spart/bpart "
- "(%zd/%zd/%zd/%zd MB)",
+ "Allocating %zd/%zd/%zd/%zd/%zd foreign part/gpart/spart/bpart/sink "
+ "(%zd/%zd/%zd/%zd/%zd MB)",
s->size_parts_foreign, s->size_gparts_foreign, s->size_sparts_foreign,
- s->size_bparts_foreign,
+ s->size_bparts_foreign, s->size_sinks_foreign,
s->size_parts_foreign * sizeof(struct part) / (1024 * 1024),
s->size_gparts_foreign * sizeof(struct gpart) / (1024 * 1024),
s->size_sparts_foreign * sizeof(struct spart) / (1024 * 1024),
- s->size_bparts_foreign * sizeof(struct bpart) / (1024 * 1024));
+ s->size_bparts_foreign * sizeof(struct bpart) / (1024 * 1024),
+ s->size_sinks_foreign * sizeof(struct sink) / (1024 * 1024));
if ((s->size_parts_foreign - old_size_parts_foreign) > 0 ||
(s->size_gparts_foreign - old_size_gparts_foreign) > 0 ||
(s->size_sparts_foreign - old_size_sparts_foreign) > 0 ||
- (s->size_bparts_foreign - old_size_bparts_foreign) > 0) {
+ (s->size_bparts_foreign - old_size_bparts_foreign) > 0 ||
+ (s->size_sinks_foreign - old_size_sinks_foreign) > 0) {
message(
- "Re-allocations %zd/%zd/%zd/%zd part/gpart/spart/bpart "
- "(%zd/%zd/%zd/%zd MB)",
+ "Re-allocations %zd/%zd/%zd/%zd/%zd part/gpart/spart/bpart/sink "
+ "(%zd/%zd/%zd/%zd/%zd MB)",
(s->size_parts_foreign - old_size_parts_foreign),
(s->size_gparts_foreign - old_size_gparts_foreign),
(s->size_sparts_foreign - old_size_sparts_foreign),
(s->size_bparts_foreign - old_size_bparts_foreign),
+ (s->size_sinks_foreign - old_size_sinks_foreign),
(s->size_parts_foreign - old_size_parts_foreign) *
sizeof(struct part) / (1024 * 1024),
(s->size_gparts_foreign - old_size_gparts_foreign) *
@@ -900,7 +927,9 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
(s->size_sparts_foreign - old_size_sparts_foreign) *
sizeof(struct spart) / (1024 * 1024),
(s->size_bparts_foreign - old_size_bparts_foreign) *
- sizeof(struct bpart) / (1024 * 1024));
+ sizeof(struct bpart) / (1024 * 1024),
+ (s->size_sinks_foreign - old_size_sinks_foreign) *
+ sizeof(struct sink) / (1024 * 1024));
}
}
@@ -909,6 +938,7 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
struct gpart *gparts = s->gparts_foreign;
struct spart *sparts = s->sparts_foreign;
struct bpart *bparts = s->bparts_foreign;
+ struct sink *sinks = s->sinks_foreign;
for (int k = 0; k < nr_proxies; k++) {
for (int j = 0; j < e->proxies[k].nr_cells_in; j++) {
@@ -940,6 +970,14 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
cell_link_bparts(e->proxies[k].cells_in[j], bparts);
bparts = &bparts[e->proxies[k].cells_in[j]->black_holes.count];
}
+
+ if (!fof && with_sinks) {
+
+ /* For sinks, we just use the numbers in the top-level cells */
+ cell_link_sinks(e->proxies[k].cells_in[j], sinks);
+ sinks =
+ &sinks[e->proxies[k].cells_in[j]->sinks.count + space_extra_sinks];
+ }
}
}
@@ -948,6 +986,7 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
s->nr_gparts_foreign = gparts - s->gparts_foreign;
s->nr_sparts_foreign = sparts - s->sparts_foreign;
s->nr_bparts_foreign = bparts - s->bparts_foreign;
+ s->nr_sinks_foreign = sinks - s->sinks_foreign;
if (e->verbose)
message("Recursively linking foreign arrays took %.3f %s.",
@@ -1135,8 +1174,15 @@ int engine_estimate_nr_tasks(const struct engine *e) {
#endif
}
if (e->policy & engine_policy_sinks) {
- /* 1 drift, 2 kicks, 1 time-step, 1 sink formation */
- n1 += 5;
+ /* 1 drift, 2 kicks, 1 time-step, 1 sink formation | 5
+ density: 1 self + 13 pairs | 14
+ swallow: 1 self + 13 pairs | 14
+ do_gas_swallow: 1 self + 13 pairs | 14
+ do_sink_swallow: 1 self + 13 pairs | 14
+ ghosts: density_ghost, sink_ghost_1, sink_ghost_2 | 3
+ implicit: sink_in, sink_out | 2 */
+ n1 += 66;
+ n2 += 3;
if (e->policy & engine_policy_stars) {
/* 1 star formation */
n1 += 1;
@@ -1205,7 +1251,8 @@ int engine_estimate_nr_tasks(const struct engine *e) {
struct cell *c = &e->s->cells_top[k];
/* Any cells with particles will have tasks (local & foreign). */
- int nparts = c->hydro.count + c->grav.count + c->stars.count;
+ int nparts = c->hydro.count + c->grav.count + c->stars.count +
+ c->black_holes.count + c->sinks.count;
if (nparts > 0) {
ntop++;
ncells++;
diff --git a/src/engine.h b/src/engine.h
index 50e98011c29ff3e2cd731ebce5fcce6b8c5df085..81c6fa07109833dcabc906e435bfb6edc1fe9a0e 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -749,7 +749,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
size_t *Ngpart, const size_t offset_sparts,
const int *ind_spart, size_t *Nspart,
const size_t offset_bparts, const int *ind_bpart,
- size_t *Nbpart);
+ size_t *Nbpart, const size_t offset_sinks,
+ const int *ind_sink, size_t *Nsink);
void engine_rebuild(struct engine *e, int redistributed, int clean_h_values);
void engine_repartition(struct engine *e);
void engine_repartition_trigger(struct engine *e);
diff --git a/src/engine_config.c b/src/engine_config.c
index 5e6c4eb98ccc427e76f771070632c082ec9683d0..74fc0d20d6dd12fc90c17d09505c092f4da148c9 100644
--- a/src/engine_config.c
+++ b/src/engine_config.c
@@ -950,6 +950,8 @@ void engine_config(int restart, int fof, struct engine *e,
params, "Scheduler:cell_extra_gparts", space_extra_gparts_default);
space_extra_bparts = parser_get_opt_param_int(
params, "Scheduler:cell_extra_bparts", space_extra_bparts_default);
+ space_extra_sinks = parser_get_opt_param_int(
+ params, "Scheduler:cell_extra_sinks", space_extra_sinks_default);
/* Do we want any spare particles for on the fly creation?
This condition should be the same than in space.c */
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index 5592e384ec9d4fc2a1609ad5744ea1b46397f8ad..e18cf26e5a1073b6f68123875d43e1f64968a908 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -365,14 +365,13 @@ void engine_addtasks_send_stars(struct engine *e, struct cell *ci,
struct cell *cj, struct task *t_density,
struct task *t_prep2, struct task *t_sf_counts,
const int with_star_formation) {
-#ifdef SWIFT_DEBUG_CHECKS
+#ifdef WITH_MPI
+#if !defined(SWIFT_DEBUG_CHECKS)
if (e->policy & engine_policy_sinks && e->policy & engine_policy_stars) {
- error("TODO");
+ error("TODO: Star formation sink over MPI");
}
#endif
-#ifdef WITH_MPI
-
struct link *l = NULL;
struct scheduler *s = &e->sched;
const int nodeID = cj->nodeID;
@@ -916,13 +915,13 @@ void engine_addtasks_recv_stars(struct engine *e, struct cell *c,
struct task *t_sf_counts,
struct task *const tend,
const int with_star_formation) {
-#ifdef SWIFT_DEBUG_CHECKS
+#ifdef WITH_MPI
+#if !defined(SWIFT_DEBUG_CHECKS)
if (e->policy & engine_policy_sinks && e->policy & engine_policy_stars) {
- error("TODO");
+ error("TODO: Star formation sink over MPI");
}
#endif
-#ifdef WITH_MPI
struct scheduler *s = &e->sched;
/* Early abort (are we below the level where tasks are)? */
@@ -4378,9 +4377,9 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
const int with_rt = (e->policy & engine_policy_rt);
struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(WITH_MPI) && !defined(SWIFT_DEBUG_CHECKS)
if (e->policy & engine_policy_sinks) {
- error("TODO");
+ error("TODO: Sink MPI tasks are not implemented yet!");
}
#endif
@@ -4451,9 +4450,9 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
const int with_rt = (e->policy & engine_policy_rt);
struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(WITH_MPI) && !defined(SWIFT_DEBUG_CHECKS)
if (e->policy & engine_policy_sinks) {
- error("TODO");
+ error("TODO: Sink MPI tasks are not implemented yet!");
}
#endif
diff --git a/src/engine_redistribute.c b/src/engine_redistribute.c
index bb8b7b6ead70796990c7ecb0311c2a65e4959bf1..cf22b146c4040d1a62c25f75b85a9975271de037 100644
--- a/src/engine_redistribute.c
+++ b/src/engine_redistribute.c
@@ -319,6 +319,14 @@ void ENGINE_REDISTRIBUTE_DEST_MAPPER(gpart);
*/
void ENGINE_REDISTRIBUTE_DEST_MAPPER(bpart);
+/**
+ * @brief Accumulate the counts of sink particles per cell.
+ * Threadpool helper for accumulating the counts of particles per cell.
+ *
+ * sink version.
+ */
+void ENGINE_REDISTRIBUTE_DEST_MAPPER(sink);
+
#endif /* redist_mapper_data */
#ifdef WITH_MPI /* savelink_mapper_data */
@@ -399,12 +407,22 @@ void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(bpart, 1);
void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(bpart, 0);
#endif
+/**
+ * @brief Save position of sink-gpart links.
+ * Threadpool helper for accumulating the counts of particles per cell.
+ */
+#ifdef SWIFT_DEBUG_CHECKS
+void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(sink, 1);
+#else
+void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(sink, 0);
+#endif
+
#endif /* savelink_mapper_data */
#ifdef WITH_MPI /* relink_mapper_data */
-/* Support for relinking parts, gparts, sparts and bparts after moving between
- * nodes. */
+/* Support for relinking parts, gparts, sparts, bparts and sinks after moving
+ * between nodes. */
struct relink_mapper_data {
int nodeID;
int nr_nodes;
@@ -412,6 +430,7 @@ struct relink_mapper_data {
int *s_counts;
int *g_counts;
int *b_counts;
+ int *sink_counts;
struct space *s;
};
@@ -435,6 +454,7 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
int *g_counts = mydata->g_counts;
int *s_counts = mydata->s_counts;
int *b_counts = mydata->b_counts;
+ int *sink_counts = mydata->sink_counts;
struct space *s = mydata->s;
for (int i = 0; i < num_elements; i++) {
@@ -446,12 +466,14 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
size_t offset_gparts = 0;
size_t offset_sparts = 0;
size_t offset_bparts = 0;
+ size_t offset_sinks = 0;
for (int n = 0; n < node; n++) {
int ind_recv = n * nr_nodes + nodeID;
offset_parts += counts[ind_recv];
offset_gparts += g_counts[ind_recv];
offset_sparts += s_counts[ind_recv];
offset_bparts += b_counts[ind_recv];
+ offset_sinks += sink_counts[ind_recv];
}
/* Number of gparts sent from this node. */
@@ -493,6 +515,17 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
s->gparts[k].id_or_neg_offset = -partner_index;
s->bparts[partner_index].gpart = &s->gparts[k];
}
+
+ /* Does this gpart have a sink partner ? */
+ else if (s->gparts[k].type == swift_type_sink) {
+
+ const ptrdiff_t partner_index =
+ offset_sinks - s->gparts[k].id_or_neg_offset;
+
+ /* Re-link */
+ s->gparts[k].id_or_neg_offset = -partner_index;
+ s->sinks[partner_index].gpart = &s->gparts[k];
+ }
}
}
}
@@ -519,13 +552,6 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
void engine_redistribute(struct engine *e) {
#ifdef WITH_MPI
-#ifdef SWIFT_DEBUG_CHECKS
- const int nr_sinks_new = 0;
-#endif
- if (e->policy & engine_policy_sinks) {
- error("Not implemented yet");
- }
-
const int nr_nodes = e->nr_nodes;
const int nodeID = e->nodeID;
struct space *s = e->s;
@@ -536,12 +562,14 @@ void engine_redistribute(struct engine *e) {
struct gpart *gparts = s->gparts;
struct spart *sparts = s->sparts;
struct bpart *bparts = s->bparts;
+ struct sink *sinks = s->sinks;
ticks tic = getticks();
size_t nr_parts = s->nr_parts;
size_t nr_gparts = s->nr_gparts;
size_t nr_sparts = s->nr_sparts;
size_t nr_bparts = s->nr_bparts;
+ size_t nr_sinks = s->nr_sinks;
/* Start by moving inhibited particles to the end of the arrays */
for (size_t k = 0; k < nr_parts; /* void */) {
@@ -609,6 +637,27 @@ void engine_redistribute(struct engine *e) {
}
}
+ /* Now move inhibited sink particles to the end of the arrays */
+ for (size_t k = 0; k < nr_sinks; /* void */) {
+ if (sinks[k].time_bin == time_bin_inhibited ||
+ sinks[k].time_bin == time_bin_not_created) {
+ nr_sinks -= 1;
+
+ /* Swap the particle */
+ memswap(&s->sinks[k], &s->sinks[nr_sinks], sizeof(struct sink));
+
+ /* Swap the link with the gpart */
+ if (s->sinks[k].gpart != NULL) {
+ s->sinks[k].gpart->id_or_neg_offset = -k;
+ }
+ if (s->sinks[nr_sinks].gpart != NULL) {
+ s->sinks[nr_sinks].gpart->id_or_neg_offset = -nr_sinks;
+ }
+ } else {
+ k++;
+ }
+ }
+
/* Finally do the same with the gravity particles */
for (size_t k = 0; k < nr_gparts; /* void */) {
if (gparts[k].time_bin == time_bin_inhibited ||
@@ -626,6 +675,8 @@ void engine_redistribute(struct engine *e) {
s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
} else if (s->gparts[k].type == swift_type_black_hole) {
s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
+ } else if (s->gparts[k].type == swift_type_sink) {
+ s->sinks[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
}
if (s->gparts[nr_gparts].type == swift_type_gas) {
@@ -637,6 +688,9 @@ void engine_redistribute(struct engine *e) {
} else if (s->gparts[nr_gparts].type == swift_type_black_hole) {
s->bparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
&s->gparts[nr_gparts];
+ } else if (s->gparts[nr_gparts].type == swift_type_sink) {
+ s->sinks[-s->gparts[nr_sinks].id_or_neg_offset].gpart =
+ &s->gparts[nr_gparts];
}
} else {
k++;
@@ -857,6 +911,73 @@ void engine_redistribute(struct engine *e) {
}
swift_free("b_dest", b_dest);
+ /* Get destination of each sink-particle */
+ int *sink_counts;
+ if ((sink_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
+ error("Failed to allocate sink_counts temporary buffer.");
+
+ int *sink_dest;
+ if ((sink_dest = (int *)swift_malloc("sink_dest", sizeof(int) * nr_sinks)) ==
+ NULL)
+ error("Failed to allocate sink_dest temporary buffer.");
+
+ redist_data.counts = sink_counts;
+ redist_data.dest = sink_dest;
+ redist_data.base = (void *)sinks;
+
+ threadpool_map(&e->threadpool, engine_redistribute_dest_mapper_sink, sinks,
+ nr_sinks, sizeof(struct sink), threadpool_auto_chunk_size,
+ &redist_data);
+
+ /* Sort the particles according to their cell index. */
+ if (nr_sinks > 0)
+ space_sinks_sort(s->sinks, sink_dest, &sink_counts[nodeID * nr_nodes],
+ nr_nodes, 0);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Verify that the sink have been sorted correctly. */
+ for (size_t k = 0; k < nr_sinks; k++) {
+ const struct sink *sink = &s->sinks[k];
+
+ if (sink->time_bin == time_bin_inhibited)
+ error("Inhibited particle found after sorting!");
+
+ if (sink->time_bin == time_bin_not_created)
+ error("Inhibited particle found after sorting!");
+
+ /* New cell index */
+ const int new_cid =
+ cell_getid(s->cdim, sink->x[0] * s->iwidth[0],
+ sink->x[1] * s->iwidth[1], sink->x[2] * s->iwidth[2]);
+
+ /* New cell of this sink */
+ const struct cell *c = &s->cells_top[new_cid];
+ const int new_node = c->nodeID;
+
+ if (sink_dest[k] != new_node)
+ error("sink's new node index not matching sorted index.");
+
+ if (sink->x[0] < c->loc[0] || sink->x[0] > c->loc[0] + c->width[0] ||
+ sink->x[1] < c->loc[1] || sink->x[1] > c->loc[1] + c->width[1] ||
+ sink->x[2] < c->loc[2] || sink->x[2] > c->loc[2] + c->width[2])
+ error("sink not sorted into the right top-level cell!");
+ }
+#endif
+
+ /* We need to re-link the gpart partners of sinks. */
+ if (nr_sinks > 0) {
+
+ struct savelink_mapper_data savelink_data;
+ savelink_data.nr_nodes = nr_nodes;
+ savelink_data.counts = sink_counts;
+ savelink_data.parts = (void *)sinks;
+ savelink_data.nodeID = nodeID;
+ threadpool_map(&e->threadpool, engine_redistribute_savelink_mapper_sink,
+ nodes, nr_nodes, sizeof(int), threadpool_auto_chunk_size,
+ &savelink_data);
+ }
+ swift_free("sink_dest", sink_dest);
+
/* Get destination of each g-particle */
int *g_counts;
if ((g_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
@@ -932,22 +1053,30 @@ void engine_redistribute(struct engine *e) {
MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS)
error("Failed to allreduce bparticle transfer counts.");
+ /* Get all the sink_counts from all the nodes. */
+ if (MPI_Allreduce(MPI_IN_PLACE, sink_counts, nr_nodes * nr_nodes, MPI_INT,
+ MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS)
+ error("Failed to allreduce sink particle transfer counts.");
+
/* Report how many particles will be moved. */
if (e->verbose) {
if (e->nodeID == 0) {
- size_t total = 0, g_total = 0, s_total = 0, b_total = 0;
- size_t unmoved = 0, g_unmoved = 0, s_unmoved = 0, b_unmoved = 0;
+ size_t total = 0, g_total = 0, s_total = 0, b_total = 0, sink_total = 0;
+ size_t unmoved = 0, g_unmoved = 0, s_unmoved = 0, b_unmoved = 0,
+ sink_unmoved = 0;
for (int p = 0, r = 0; p < nr_nodes; p++) {
for (int n = 0; n < nr_nodes; n++) {
total += counts[r];
g_total += g_counts[r];
s_total += s_counts[r];
b_total += b_counts[r];
+ sink_total += sink_counts[r];
if (p == n) {
unmoved += counts[r];
g_unmoved += g_counts[r];
s_unmoved += s_counts[r];
b_unmoved += b_counts[r];
+ sink_unmoved += sink_counts[r];
}
r++;
}
@@ -967,14 +1096,19 @@ void engine_redistribute(struct engine *e) {
message("%ld of %ld (%.2f%%) of b-particles moved", b_total - b_unmoved,
b_total,
100.0 * (double)(b_total - b_unmoved) / (double)b_total);
+ if (sink_total > 0)
+ message(
+ "%ld of %ld (%.2f%%) of sink-particles moved",
+ sink_total - sink_unmoved, sink_total,
+ 100.0 * (double)(sink_total - sink_unmoved) / (double)sink_total);
}
}
- /* Now each node knows how many parts, sparts, bparts, and gparts will be
- * transferred to every other node. Get the new numbers of particles for this
- * node. */
+ /* Now each node knows how many parts, sparts, bparts, sinks and gparts will
+ * be transferred to every other node. Get the new numbers of particles for
+ * this node. */
size_t nr_parts_new = 0, nr_gparts_new = 0, nr_sparts_new = 0,
- nr_bparts_new = 0;
+ nr_bparts_new = 0, nr_sinks_new = 0;
for (int k = 0; k < nr_nodes; k++)
nr_parts_new += counts[k * nr_nodes + nodeID];
for (int k = 0; k < nr_nodes; k++)
@@ -983,6 +1117,8 @@ void engine_redistribute(struct engine *e) {
nr_sparts_new += s_counts[k * nr_nodes + nodeID];
for (int k = 0; k < nr_nodes; k++)
nr_bparts_new += b_counts[k * nr_nodes + nodeID];
+ for (int k = 0; k < nr_nodes; k++)
+ nr_sinks_new += sink_counts[k * nr_nodes + nodeID];
#ifdef WITH_CSDS
const int initial_redistribute = e->ti_current == 0;
@@ -992,6 +1128,7 @@ void engine_redistribute(struct engine *e) {
size_t spart_offset = 0;
size_t gpart_offset = 0;
size_t bpart_offset = 0;
+ size_t sink_offset = 0;
for (int i = 0; i < nr_nodes; i++) {
const size_t c_ind = engine_rank * nr_nodes + i;
@@ -1002,6 +1139,7 @@ void engine_redistribute(struct engine *e) {
spart_offset += s_counts[c_ind];
gpart_offset += g_counts[c_ind];
bpart_offset += b_counts[c_ind];
+ sink_offset += sink_counts[c_ind];
continue;
}
@@ -1023,11 +1161,17 @@ void engine_redistribute(struct engine *e) {
error("TODO");
}
+ /* Log the sinks */
+ if (sink_counts[c_ind] > 0) {
+ error("TODO");
+ }
+
/* Update the counters */
part_offset += counts[c_ind];
spart_offset += s_counts[c_ind];
gpart_offset += g_counts[c_ind];
bpart_offset += b_counts[c_ind];
+ sink_offset += sink_counts[c_ind];
}
}
#endif
@@ -1081,6 +1225,15 @@ void engine_redistribute(struct engine *e) {
s->nr_bparts = nr_bparts_new;
s->size_bparts = engine_redistribute_alloc_margin * nr_bparts_new;
+ /* Sink particles. */
+ new_parts = engine_do_redistribute(
+ "sinks", sink_counts, (char *)s->sinks, nr_sinks_new, sizeof(struct sink),
+ sink_align, sink_mpi_type, nr_nodes, nodeID, e->syncredist);
+ swift_free("sinks", s->sinks);
+ s->sinks = (struct sink *)new_parts;
+ s->nr_sinks = nr_sinks_new;
+ s->size_sinks = engine_redistribute_alloc_margin * nr_sinks_new;
+
/* All particles have now arrived. Time for some final operations on the
stuff we just received */
@@ -1090,6 +1243,7 @@ void engine_redistribute(struct engine *e) {
size_t spart_offset = 0;
size_t gpart_offset = 0;
size_t bpart_offset = 0;
+ size_t sink_offset = 0;
for (int i = 0; i < nr_nodes; i++) {
const size_t c_ind = i * nr_nodes + engine_rank;
@@ -1100,6 +1254,7 @@ void engine_redistribute(struct engine *e) {
spart_offset += s_counts[c_ind];
gpart_offset += g_counts[c_ind];
bpart_offset += b_counts[c_ind];
+ sink_offset += sink_counts[c_ind];
continue;
}
@@ -1121,11 +1276,17 @@ void engine_redistribute(struct engine *e) {
error("TODO");
}
+ /* Log the sinks */
+ if (sink_counts[c_ind] > 0) {
+ error("TODO");
+ }
+
/* Update the counters */
part_offset += counts[c_ind];
spart_offset += s_counts[c_ind];
gpart_offset += g_counts[c_ind];
bpart_offset += b_counts[c_ind];
+ sink_offset += sink_counts[c_ind];
}
}
#endif
@@ -1139,6 +1300,7 @@ void engine_redistribute(struct engine *e) {
relink_data.g_counts = g_counts;
relink_data.s_counts = s_counts;
relink_data.b_counts = b_counts;
+ relink_data.sink_counts = sink_counts;
relink_data.nodeID = nodeID;
relink_data.nr_nodes = nr_nodes;
@@ -1151,6 +1313,7 @@ void engine_redistribute(struct engine *e) {
free(g_counts);
free(s_counts);
free(b_counts);
+ free(sink_counts);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that all parts are in the right place. */
@@ -1186,6 +1349,15 @@ void engine_redistribute(struct engine *e) {
error("Received b-particle (%zu) that does not belong here (nodeID=%i).",
k, cells[cid].nodeID);
}
+ for (size_t k = 0; k < nr_sinks_new; k++) {
+ const int cid = cell_getid(s->cdim, s->sinks[k].x[0] * s->iwidth[0],
+ s->sinks[k].x[1] * s->iwidth[1],
+ s->sinks[k].x[2] * s->iwidth[2]);
+ if (cells[cid].nodeID != nodeID)
+ error(
+ "Received sink-particle (%zu) that does not belong here (nodeID=%i).",
+ k, cells[cid].nodeID);
+ }
/* Verify that the links are correct */
part_verify_links(s->parts, s->gparts, s->sinks, s->sparts, s->bparts,
@@ -1200,10 +1372,11 @@ void engine_redistribute(struct engine *e) {
for (int k = 0; k < nr_cells; k++)
if (cells[k].nodeID == nodeID) my_cells += 1;
message(
- "node %i now has %zu parts, %zu sparts, %zu bparts and %zu gparts in "
+ "node %i now has %zu parts, %zu sparts, %zu bparts, %zu sinks and %zu "
+ "gparts in "
"%i cells.",
- nodeID, nr_parts_new, nr_sparts_new, nr_bparts_new, nr_gparts_new,
- my_cells);
+ nodeID, nr_parts_new, nr_sparts_new, nr_bparts_new, nr_sinks_new,
+ nr_gparts_new, my_cells);
}
/* Flag that we do not have any extra particles any more */
@@ -1211,6 +1384,7 @@ void engine_redistribute(struct engine *e) {
s->nr_extra_gparts = 0;
s->nr_extra_sparts = 0;
s->nr_extra_bparts = 0;
+ s->nr_extra_sinks = 0;
/* Flag that a redistribute has taken place */
e->step_props |= engine_step_prop_redistribute;
diff --git a/src/engine_strays.c b/src/engine_strays.c
index d55909c73e08c5b594d6ab30662f3b6e64d92a25..2ecd2576fcb0091534aaf64e83550eef9d82234e 100644
--- a/src/engine_strays.c
+++ b/src/engine_strays.c
@@ -33,6 +33,13 @@
/* Local headers. */
#include "proxy.h"
+#ifdef WITH_MPI
+/* Number of particle types to wait for after launching the proxies. We have
+ parts, xparts, gparts, sparts, bparts and sinks to exchange, hence 6 types.
+ */
+#define MPI_REQUEST_NUMBER_PARTICLE_TYPES 6
+#endif
+
/**
* @brief Exchange straying particles with other nodes.
*
@@ -57,18 +64,22 @@
* @param ind_bpart The foreign #cell ID of each bpart.
* @param Nbpart The number of stray bparts, contains the number of bparts
* received on return.
+ * @param offset_sinks The index in the sinks array as of which the foreign
+ * parts reside (i.e. the current number of local #sink).
+ * @param ind_sink The foreign #cell ID of each sink.
+ * @param Nsink The number of stray sinks, contains the number of sinks
+ * received on return.
*
* Note that this function does not mess-up the linkage between parts and
* gparts, i.e. the received particles have correct linkeage.
*/
-void engine_exchange_strays(struct engine *e, const size_t offset_parts,
- const int *restrict ind_part, size_t *Npart,
- const size_t offset_gparts,
- const int *restrict ind_gpart, size_t *Ngpart,
- const size_t offset_sparts,
- const int *restrict ind_spart, size_t *Nspart,
- const size_t offset_bparts,
- const int *restrict ind_bpart, size_t *Nbpart) {
+void engine_exchange_strays(
+ struct engine *e, const size_t offset_parts, const int *restrict ind_part,
+ size_t *Npart, const size_t offset_gparts, const int *restrict ind_gpart,
+ size_t *Ngpart, const size_t offset_sparts, const int *restrict ind_spart,
+ size_t *Nspart, const size_t offset_bparts, const int *restrict ind_bpart,
+ size_t *Nbpart, const size_t offset_sinks, const int *restrict ind_sink,
+ size_t *Nsink) {
#ifdef WITH_MPI
struct space *s = e->s;
@@ -80,6 +91,7 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
e->proxies[k].nr_gparts_out = 0;
e->proxies[k].nr_sparts_out = 0;
e->proxies[k].nr_bparts_out = 0;
+ e->proxies[k].nr_sinks_out = 0;
}
/* Put the parts into the corresponding proxies. */
@@ -204,6 +216,46 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
/* Load the bpart into the proxy */
proxy_bparts_load(&e->proxies[pid], &s->bparts[offset_bparts + k], 1);
+#ifdef WITH_CSDS
+ if (e->policy & engine_policy_csds) {
+ error("Not yet implemented.");
+ }
+#endif
+ }
+
+ /* Put the sinks into the corresponding proxies. */
+ for (size_t k = 0; k < *Nsink; k++) {
+ /* Ignore the particles we want to get rid of (inhibited, ...). */
+ if (ind_sink[k] == -1) continue;
+
+ /* Get the target node and proxy ID. */
+ const int node_id = e->s->cells_top[ind_sink[k]].nodeID;
+ if (node_id < 0 || node_id >= e->nr_nodes)
+ error("Bad node ID %i.", node_id);
+ const int pid = e->proxy_ind[node_id];
+ if (pid < 0) {
+ error(
+ "Do not have a proxy for the requested nodeID %i for part with "
+ "id=%lld, x=[%e,%e,%e].",
+ node_id, s->sinks[offset_sinks + k].id,
+ s->sinks[offset_sinks + k].x[0], s->sinks[offset_sinks + k].x[1],
+ s->sinks[offset_sinks + k].x[2]);
+ }
+
+ /* Re-link the associated gpart with the buffer offset of the sink. */
+ if (s->sinks[offset_sinks + k].gpart != NULL) {
+ s->sinks[offset_sinks + k].gpart->id_or_neg_offset =
+ -e->proxies[pid].nr_sinks_out;
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (s->sinks[offset_sinks + k].time_bin == time_bin_inhibited)
+ error("Attempting to exchange an inhibited particle");
+#endif
+
+ /* Load the sink into the proxy */
+ proxy_sinks_load(&e->proxies[pid], &s->sinks[offset_sinks + k], 1);
+
#ifdef WITH_CSDS
if (e->policy & engine_policy_csds) {
error("Not yet implemented.");
@@ -252,8 +304,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
}
/* Launch the proxies. */
- MPI_Request reqs_in[5 * engine_maxproxies];
- MPI_Request reqs_out[5 * engine_maxproxies];
+ MPI_Request reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * engine_maxproxies];
+ MPI_Request reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * engine_maxproxies];
for (int k = 0; k < e->nr_proxies; k++) {
proxy_parts_exchange_first(&e->proxies[k]);
reqs_in[k] = e->proxies[k].req_parts_count_in;
@@ -281,18 +333,23 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
int count_gparts_in = 0;
int count_sparts_in = 0;
int count_bparts_in = 0;
+ int count_sinks_in = 0;
for (int k = 0; k < e->nr_proxies; k++) {
count_parts_in += e->proxies[k].nr_parts_in;
count_gparts_in += e->proxies[k].nr_gparts_in;
count_sparts_in += e->proxies[k].nr_sparts_in;
count_bparts_in += e->proxies[k].nr_bparts_in;
+ count_sinks_in += e->proxies[k].nr_sinks_in;
+ message("Counting entering particles, k = %i, nr_proxies = %d", k,
+ e->nr_proxies);
}
if (e->verbose) {
message(
- "sent out %zu/%zu/%zu/%zu parts/gparts/sparts/bparts, got %i/%i/%i/%i "
+ "sent out %zu/%zu/%zu/%zu/%zu parts/gparts/sparts/bparts/sinks, got "
+ "%i/%i/%i/%i/%i "
"back.",
- *Npart, *Ngpart, *Nspart, *Nbpart, count_parts_in, count_gparts_in,
- count_sparts_in, count_bparts_in);
+ *Npart, *Ngpart, *Nspart, *Nbpart, *Nsink, count_parts_in,
+ count_gparts_in, count_sparts_in, count_bparts_in, count_sinks_in);
}
/* Reallocate the particle arrays if necessary */
@@ -356,6 +413,24 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
}
}
+ if (offset_sinks + count_sinks_in > s->size_sinks) {
+ s->size_sinks = (offset_sinks + count_sinks_in) * engine_parts_size_grow;
+ struct sink *sinks_new = NULL;
+ if (swift_memalign("sinks", (void **)&sinks_new, sink_align,
+ sizeof(struct sink) * s->size_sinks) != 0)
+ error("Failed to allocate new sink data.");
+ memcpy(sinks_new, s->sinks, sizeof(struct sink) * offset_sinks);
+ swift_free("sinks", s->sinks);
+ s->sinks = sinks_new;
+
+ /* Reset the links */
+ for (size_t k = 0; k < offset_sinks; k++) {
+ if (s->sinks[k].gpart != NULL) {
+ s->sinks[k].gpart->id_or_neg_offset = -k;
+ }
+ }
+ }
+
if (offset_gparts + count_gparts_in > s->size_gparts) {
s->size_gparts = (offset_gparts + count_gparts_in) * engine_parts_size_grow;
struct gpart *gparts_new = NULL;
@@ -374,6 +449,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
} else if (s->gparts[k].type == swift_type_black_hole) {
s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
+ } else if (s->gparts[k].type == swift_type_sink) {
+ s->sinks[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
}
}
}
@@ -382,82 +459,113 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
int nr_in = 0, nr_out = 0;
for (int k = 0; k < e->nr_proxies; k++) {
if (e->proxies[k].nr_parts_in > 0) {
- reqs_in[5 * k] = e->proxies[k].req_parts_in;
- reqs_in[5 * k + 1] = e->proxies[k].req_xparts_in;
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
+ e->proxies[k].req_parts_in;
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] =
+ e->proxies[k].req_xparts_in;
nr_in += 2;
} else {
- reqs_in[5 * k] = reqs_in[5 * k + 1] = MPI_REQUEST_NULL;
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_gparts_in > 0) {
- reqs_in[5 * k + 2] = e->proxies[k].req_gparts_in;
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] =
+ e->proxies[k].req_gparts_in;
nr_in += 1;
} else {
- reqs_in[5 * k + 2] = MPI_REQUEST_NULL;
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_sparts_in > 0) {
- reqs_in[5 * k + 3] = e->proxies[k].req_sparts_in;
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] =
+ e->proxies[k].req_sparts_in;
nr_in += 1;
} else {
- reqs_in[5 * k + 3] = MPI_REQUEST_NULL;
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_bparts_in > 0) {
- reqs_in[5 * k + 4] = e->proxies[k].req_bparts_in;
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] =
+ e->proxies[k].req_bparts_in;
+ nr_in += 1;
+ } else {
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] = MPI_REQUEST_NULL;
+ }
+ if (e->proxies[k].nr_sinks_in > 0) {
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] =
+ e->proxies[k].req_sinks_in;
nr_in += 1;
} else {
- reqs_in[5 * k + 4] = MPI_REQUEST_NULL;
+ reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_parts_out > 0) {
- reqs_out[5 * k] = e->proxies[k].req_parts_out;
- reqs_out[5 * k + 1] = e->proxies[k].req_xparts_out;
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
+ e->proxies[k].req_parts_out;
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] =
+ e->proxies[k].req_xparts_out;
nr_out += 2;
} else {
- reqs_out[5 * k] = reqs_out[5 * k + 1] = MPI_REQUEST_NULL;
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] =
+ MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_gparts_out > 0) {
- reqs_out[5 * k + 2] = e->proxies[k].req_gparts_out;
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] =
+ e->proxies[k].req_gparts_out;
nr_out += 1;
} else {
- reqs_out[5 * k + 2] = MPI_REQUEST_NULL;
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_sparts_out > 0) {
- reqs_out[5 * k + 3] = e->proxies[k].req_sparts_out;
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] =
+ e->proxies[k].req_sparts_out;
nr_out += 1;
} else {
- reqs_out[5 * k + 3] = MPI_REQUEST_NULL;
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_bparts_out > 0) {
- reqs_out[5 * k + 4] = e->proxies[k].req_bparts_out;
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] =
+ e->proxies[k].req_bparts_out;
nr_out += 1;
} else {
- reqs_out[5 * k + 4] = MPI_REQUEST_NULL;
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] = MPI_REQUEST_NULL;
+ }
+ if (e->proxies[k].nr_sinks_out > 0) {
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] =
+ e->proxies[k].req_sinks_out;
+ nr_out += 1;
+ } else {
+ reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] = MPI_REQUEST_NULL;
}
}
/* Wait for each part array to come in and collect the new
parts from the proxies. */
- int count_parts = 0, count_gparts = 0, count_sparts = 0, count_bparts = 0;
+ int count_parts = 0, count_gparts = 0, count_sparts = 0, count_bparts = 0,
+ count_sinks = 0;
for (int k = 0; k < nr_in; k++) {
int err, pid;
- if ((err = MPI_Waitany(5 * e->nr_proxies, reqs_in, &pid,
- MPI_STATUS_IGNORE)) != MPI_SUCCESS) {
+ if ((err = MPI_Waitany(MPI_REQUEST_NUMBER_PARTICLE_TYPES * e->nr_proxies,
+ reqs_in, &pid, MPI_STATUS_IGNORE)) != MPI_SUCCESS) {
char buff[MPI_MAX_ERROR_STRING];
int res;
MPI_Error_string(err, buff, &res);
error("MPI_Waitany failed (%s).", buff);
}
if (pid == MPI_UNDEFINED) break;
- // message( "request from proxy %i has arrived." , pid / 5 );
- pid = 5 * (pid / 5);
+ // message( "request from proxy %i has arrived." , pid /
+ // MPI_REQUEST_NUMBER_PARTICLE_TYPES );
+ pid = MPI_REQUEST_NUMBER_PARTICLE_TYPES *
+ (pid / MPI_REQUEST_NUMBER_PARTICLE_TYPES);
/* If all the requests for a given proxy have arrived... */
if (reqs_in[pid + 0] == MPI_REQUEST_NULL &&
reqs_in[pid + 1] == MPI_REQUEST_NULL &&
reqs_in[pid + 2] == MPI_REQUEST_NULL &&
reqs_in[pid + 3] == MPI_REQUEST_NULL &&
- reqs_in[pid + 4] == MPI_REQUEST_NULL) {
+ reqs_in[pid + 4] == MPI_REQUEST_NULL &&
+ reqs_in[pid + 5] == MPI_REQUEST_NULL) {
/* Copy the particle data to the part/xpart/gpart arrays. */
- struct proxy *prox = &e->proxies[pid / 5];
+ struct proxy *prox = &e->proxies[pid / MPI_REQUEST_NUMBER_PARTICLE_TYPES];
memcpy(&s->parts[offset_parts + count_parts], prox->parts_in,
sizeof(struct part) * prox->nr_parts_in);
memcpy(&s->xparts[offset_parts + count_parts], prox->xparts_in,
@@ -468,6 +576,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
sizeof(struct spart) * prox->nr_sparts_in);
memcpy(&s->bparts[offset_bparts + count_bparts], prox->bparts_in,
sizeof(struct bpart) * prox->nr_bparts_in);
+ memcpy(&s->sinks[offset_sinks + count_sinks], prox->sinks_in,
+ sizeof(struct sink) * prox->nr_sinks_in);
#ifdef WITH_CSDS
if (e->policy & engine_policy_csds) {
@@ -495,6 +605,12 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
if (prox->nr_bparts_in > 0) {
error("TODO");
}
+
+ /* Log the sinks */
+ if (prox->nr_sinks_in > 0) {
+ /* Not implemented yet */
+ error("TODO");
+ }
}
#endif
/* for (int k = offset; k < offset + count; k++)
@@ -522,6 +638,11 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
&s->bparts[offset_bparts + count_bparts - gp->id_or_neg_offset];
gp->id_or_neg_offset = s->bparts - bp;
bp->gpart = gp;
+ } else if (gp->type == swift_type_sink) {
+ struct sink *sink =
+ &s->sinks[offset_sinks + count_sinks - gp->id_or_neg_offset];
+ gp->id_or_neg_offset = s->sinks - sink;
+ sink->gpart = gp;
}
}
@@ -530,13 +651,14 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
count_gparts += prox->nr_gparts_in;
count_sparts += prox->nr_sparts_in;
count_bparts += prox->nr_bparts_in;
+ count_sinks += prox->nr_sinks_in;
}
}
/* Wait for all the sends to have finished too. */
if (nr_out > 0)
- if (MPI_Waitall(5 * e->nr_proxies, reqs_out, MPI_STATUSES_IGNORE) !=
- MPI_SUCCESS)
+ if (MPI_Waitall(MPI_REQUEST_NUMBER_PARTICLE_TYPES * e->nr_proxies, reqs_out,
+ MPI_STATUSES_IGNORE) != MPI_SUCCESS)
error("MPI_Waitall on sends failed.");
/* Free the proxy memory */
@@ -553,6 +675,7 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
*Ngpart = count_gparts;
*Nspart = count_sparts;
*Nbpart = count_bparts;
+ *Nsink = count_sinks;
#else
error("SWIFT was not compiled with MPI support.");
diff --git a/src/feedback/GEAR/feedback.c b/src/feedback/GEAR/feedback.c
index 46dfc6dea1ced3d51dc7eea992f1ea3c3657ad98..204477d474612d95419117d6d745844454c4330a 100644
--- a/src/feedback/GEAR/feedback.c
+++ b/src/feedback/GEAR/feedback.c
@@ -21,6 +21,7 @@
#include "feedback.h"
/* Local includes */
+#include "cooling.h"
#include "cosmology.h"
#include "engine.h"
#include "error.h"
@@ -49,7 +50,7 @@ void feedback_update_part(struct part* p, struct xpart* xp,
const struct pressure_floor_props* pressure_floor = e->pressure_floor_props;
/* Turn off the cooling */
- xp->cooling_data.time_last_event = e->time;
+ cooling_set_part_time_cooling_off(p, xp, e->time);
/* Update mass */
const float old_mass = hydro_get_mass(p);
diff --git a/src/part.c b/src/part.c
index 9d9c38f577f8675d7cb9367a2227bce689b892f0..f513153d61ba9b277a0b6b88ae9e4cf61c766e10 100644
--- a/src/part.c
+++ b/src/part.c
@@ -544,6 +544,7 @@ MPI_Datatype xpart_mpi_type;
MPI_Datatype gpart_mpi_type;
MPI_Datatype spart_mpi_type;
MPI_Datatype bpart_mpi_type;
+MPI_Datatype sink_mpi_type;
/**
* @brief Registers MPI particle types.
@@ -581,6 +582,11 @@ void part_create_mpi_types(void) {
MPI_Type_commit(&bpart_mpi_type) != MPI_SUCCESS) {
error("Failed to create MPI type for bparts.");
}
+ if (MPI_Type_contiguous(sizeof(struct sink) / sizeof(unsigned char), MPI_BYTE,
+ &sink_mpi_type) != MPI_SUCCESS ||
+ MPI_Type_commit(&sink_mpi_type) != MPI_SUCCESS) {
+ error("Failed to create MPI type for sink.");
+ }
}
void part_free_mpi_types(void) {
@@ -590,5 +596,6 @@ void part_free_mpi_types(void) {
MPI_Type_free(&gpart_mpi_type);
MPI_Type_free(&spart_mpi_type);
MPI_Type_free(&bpart_mpi_type);
+ MPI_Type_free(&sink_mpi_type);
}
#endif
diff --git a/src/part.h b/src/part.h
index b8fc625885a9f329633bc39476ecb7aac2228b79..c0e87bfde7c6525b7cdc05974d56fa4bb5be5bd9 100644
--- a/src/part.h
+++ b/src/part.h
@@ -171,6 +171,7 @@ extern MPI_Datatype xpart_mpi_type;
extern MPI_Datatype gpart_mpi_type;
extern MPI_Datatype spart_mpi_type;
extern MPI_Datatype bpart_mpi_type;
+extern MPI_Datatype sink_mpi_type;
void part_create_mpi_types(void);
void part_free_mpi_types(void);
diff --git a/src/proxy.c b/src/proxy.c
index 22663d42a70c75b3966d0272d93a294480057b0b..3830162ada2135d28c68b3763661a18e17bcbab0 100644
--- a/src/proxy.c
+++ b/src/proxy.c
@@ -705,12 +705,24 @@ void proxy_parts_exchange_first(struct proxy *p) {
p->buff_out[1] = p->nr_gparts_out;
p->buff_out[2] = p->nr_sparts_out;
p->buff_out[3] = p->nr_bparts_out;
- if (MPI_Isend(p->buff_out, 4, MPI_INT, p->nodeID,
- p->mynodeID * proxy_tag_shift + proxy_tag_count, MPI_COMM_WORLD,
- &p->req_parts_count_out) != MPI_SUCCESS)
+ p->buff_out[4] = p->nr_sinks_out;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ message("Number of particles out [%i , %i, %i, %i, %i]", p->nr_parts_out,
+ p->nr_gparts_out, p->nr_sparts_out, p->nr_bparts_out,
+ p->nr_sinks_out);
+#endif /* SWIFT_DEBUG_CHECKS */
+
+ if (MPI_Isend(p->buff_out, PROXY_EXCHANGE_NUMBER_PARTICLE_TYPES, MPI_INT,
+ p->nodeID, p->mynodeID * proxy_tag_shift + proxy_tag_count,
+ MPI_COMM_WORLD, &p->req_parts_count_out) != MPI_SUCCESS)
error("Failed to isend nr of parts.");
- /* message( "isent particle counts [%i, %i] from node %i to node %i." ,
- p->buff_out[0], p->buff_out[1], p->mynodeID , p->nodeID ); fflush(stdout); */
+#ifdef SWIFT_DEBUG_CHECKS
+ message("isent particle counts [%i, %i, %i, %i, %i] from node %i to node %i.",
+ p->buff_out[0], p->buff_out[1], p->buff_out[2], p->buff_out[3],
+ p->buff_out[4], p->mynodeID, p->nodeID);
+ fflush(stdout);
+#endif /* SWIFT_DEBUG_CHECKS */
/* Send the particle buffers. */
if (p->nr_parts_out > 0) {
@@ -721,20 +733,28 @@ void proxy_parts_exchange_first(struct proxy *p) {
p->mynodeID * proxy_tag_shift + proxy_tag_xparts,
MPI_COMM_WORLD, &p->req_xparts_out) != MPI_SUCCESS)
error("Failed to isend part data.");
- // message( "isent particle data (%i) to node %i." , p->nr_parts_out ,
- // p->nodeID ); fflush(stdout);
- /*for (int k = 0; k < p->nr_parts_out; k++)
+#ifdef SWIFT_DEBUG_CHECKS
+ message("isent particle data (%i) to node %i.", p->nr_parts_out, p->nodeID);
+ fflush(stdout);
+ for (int k = 0; k < p->nr_parts_out; k++)
message("sending particle %lli, x=[%.3e %.3e %.3e], h=%.3e, to node %i.",
p->parts_out[k].id, p->parts_out[k].x[0], p->parts_out[k].x[1],
- p->parts_out[k].x[2], p->parts_out[k].h, p->nodeID);*/
+ p->parts_out[k].x[2], p->parts_out[k].h, p->nodeID);
+#endif /* SWIFT_DEBUG_CHECKS */
}
if (p->nr_gparts_out > 0) {
if (MPI_Isend(p->gparts_out, p->nr_gparts_out, gpart_mpi_type, p->nodeID,
p->mynodeID * proxy_tag_shift + proxy_tag_gparts,
MPI_COMM_WORLD, &p->req_gparts_out) != MPI_SUCCESS)
error("Failed to isend gpart data.");
- // message( "isent gpart data (%i) to node %i." , p->nr_gparts_out ,
- // p->nodeID ); fflush(stdout);
+#ifdef SWIFT_DEBUG_CHECKS
+ message("isent gpart data (%i) to node %i.", p->nr_gparts_out, p->nodeID);
+ fflush(stdout);
+ for (int k = 0; k < p->nr_parts_out; k++)
+ message("sending gpart %lli, x=[%.3e %.3e %.3e], to node %i.",
+ p->gparts_out[k].id_or_neg_offset, p->gparts_out[k].x[0],
+ p->gparts_out[k].x[1], p->gparts_out[k].x[2], p->nodeID);
+#endif /* SWIFT_DEBUG_CHECKS */
}
if (p->nr_sparts_out > 0) {
@@ -742,23 +762,56 @@ void proxy_parts_exchange_first(struct proxy *p) {
p->mynodeID * proxy_tag_shift + proxy_tag_sparts,
MPI_COMM_WORLD, &p->req_sparts_out) != MPI_SUCCESS)
error("Failed to isend spart data.");
- // message( "isent spart data (%i) to node %i." , p->nr_sparts_out ,
- // p->nodeID ); fflush(stdout);
+#ifdef SWIFT_DEBUG_CHECKS
+ message("isent spart data (%i) to node %i.", p->nr_sparts_out, p->nodeID);
+ fflush(stdout);
+ for (int k = 0; k < p->nr_sparts_out; k++)
+ message("sending spart %lli, x=[%.3e %.3e %.3e], h=%.3e, to node %i.",
+ p->sparts_out[k].id, p->sparts_out[k].x[0], p->sparts_out[k].x[1],
+ p->sparts_out[k].x[2], p->sparts_out[k].h, p->nodeID);
+#endif /* SWIFT_DEBUG_CHECKS */
}
if (p->nr_bparts_out > 0) {
if (MPI_Isend(p->bparts_out, p->nr_bparts_out, bpart_mpi_type, p->nodeID,
p->mynodeID * proxy_tag_shift + proxy_tag_bparts,
MPI_COMM_WORLD, &p->req_bparts_out) != MPI_SUCCESS)
error("Failed to isend bpart data.");
- // message( "isent bpart data (%i) to node %i." , p->nr_bparts_out ,
- // p->nodeID ); fflush(stdout);
+#ifdef SWIFT_DEBUG_CHECKS
+ message("isent bpart data (%i) to node %i.", p->nr_bparts_out, p->nodeID);
+ fflush(stdout);
+ for (int k = 0; k < p->nr_bparts_out; k++)
+ message("sending bpart %lli, x=[%.3e %.3e %.3e], h=%.3e, to node %i.",
+ p->bparts_out[k].id, p->bparts_out[k].x[0], p->bparts_out[k].x[1],
+ p->bparts_out[k].x[2], p->bparts_out[k].h, p->nodeID);
+#endif /* SWIFT_DEBUG_CHECKS */
+ }
+ if (p->nr_sinks_out > 0) {
+ if (MPI_Isend(p->sinks_out, p->nr_sinks_out, sink_mpi_type, p->nodeID,
+ p->mynodeID * proxy_tag_shift + proxy_tag_sinks,
+ MPI_COMM_WORLD, &p->req_sinks_out) != MPI_SUCCESS)
+ error("Failed to isend sink data.");
+#ifdef SWIFT_DEBUG_CHECKS
+ message("isent sink data (%i) to node %i.", p->nr_sinks_out, p->nodeID);
+ fflush(stdout);
+ for (int k = 0; k < p->nr_sinks_out; k++)
+ message("sending sinks %lli, x=[%.3e %.3e %.3e], h=%.3e, to node %i.",
+ p->sinks_out[k].id, p->sinks_out[k].x[0], p->sinks_out[k].x[1],
+ p->sinks_out[k].x[2], p->sinks_out[k].r_cut, p->nodeID);
+#endif /* SWIFT_DEBUG_CHECKS */
}
/* Receive the number of particles. */
- if (MPI_Irecv(p->buff_in, 4, MPI_INT, p->nodeID,
- p->nodeID * proxy_tag_shift + proxy_tag_count, MPI_COMM_WORLD,
- &p->req_parts_count_in) != MPI_SUCCESS)
+ if (MPI_Irecv(p->buff_in, PROXY_EXCHANGE_NUMBER_PARTICLE_TYPES, MPI_INT,
+ p->nodeID, p->nodeID * proxy_tag_shift + proxy_tag_count,
+ MPI_COMM_WORLD, &p->req_parts_count_in) != MPI_SUCCESS)
error("Failed to irecv nr of parts.");
+#ifdef SWIFT_DEBUG_CHECKS
+ message(
+ "irecv particle counts [%i, %i, %i, %i, %i] from node %i, I am node %i.",
+ p->buff_in[0], p->buff_in[1], p->buff_in[2], p->buff_in[3], p->buff_in[4],
+ p->nodeID, p->mynodeID);
+ fflush(stdout);
+#endif /* SWIFT_DEBUG_CHECKS */
#else
error("SWIFT was not compiled with MPI support.");
@@ -774,6 +827,7 @@ void proxy_parts_exchange_second(struct proxy *p) {
p->nr_gparts_in = p->buff_in[1];
p->nr_sparts_in = p->buff_in[2];
p->nr_bparts_in = p->buff_in[3];
+ p->nr_sinks_in = p->buff_in[4];
/* Is there enough space in the buffers? */
if (p->nr_parts_in > p->size_parts_in) {
@@ -815,6 +869,15 @@ void proxy_parts_exchange_second(struct proxy *p) {
"bparts_in", sizeof(struct bpart) * p->size_bparts_in)) == NULL)
error("Failed to re-allocate bparts_in buffers.");
}
+ if (p->nr_sinks_in > p->size_sinks_in) {
+ do {
+ p->size_sinks_in *= proxy_buffgrow;
+ } while (p->nr_sinks_in > p->size_sinks_in);
+ swift_free("sinks_in", p->sinks_in);
+ if ((p->sinks_in = (struct sink *)swift_malloc(
+ "sinks_in", sizeof(struct sink) * p->size_sinks_in)) == NULL)
+ error("Failed to re-allocate sinks_in buffers.");
+ }
/* Receive the particle buffers. */
if (p->nr_parts_in > 0) {
@@ -825,32 +888,73 @@ void proxy_parts_exchange_second(struct proxy *p) {
p->nodeID * proxy_tag_shift + proxy_tag_xparts,
MPI_COMM_WORLD, &p->req_xparts_in) != MPI_SUCCESS)
error("Failed to irecv part data.");
- // message( "irecv particle data (%i) from node %i." , p->nr_parts_in ,
- // p->nodeID ); fflush(stdout);
+#ifdef SWIFT_DEBUG_CHECKS
+ message("irecv particle data (%i) from node %i.", p->nr_parts_in,
+ p->nodeID);
+ fflush(stdout);
+ for (int k = 0; k < p->nr_parts_in; k++)
+ message("receiving parts %lli, x=[%.3e %.3e %.3e], h=%.3e, from node %i.",
+ p->parts_in[k].id, p->parts_in[k].x[0], p->parts_in[k].x[1],
+ p->parts_in[k].x[2], p->parts_in[k].h, p->nodeID);
+#endif /* SWIFT_DEBUG_CHECKS */
}
if (p->nr_gparts_in > 0) {
if (MPI_Irecv(p->gparts_in, p->nr_gparts_in, gpart_mpi_type, p->nodeID,
p->nodeID * proxy_tag_shift + proxy_tag_gparts,
MPI_COMM_WORLD, &p->req_gparts_in) != MPI_SUCCESS)
error("Failed to irecv gpart data.");
- // message( "irecv gpart data (%i) from node %i." , p->nr_gparts_in ,
- // p->nodeID ); fflush(stdout);
+#ifdef SWIFT_DEBUG_CHECKS
+ message("irecv gpart data (%i) from node %i.", p->nr_gparts_in, p->nodeID);
+ fflush(stdout);
+ for (int k = 0; k < p->nr_gparts_in; k++)
+ message("receiving gparts %lli, x=[%.3e %.3e %.3e], from node %i.",
+ p->gparts_in[k].id_or_neg_offset, p->gparts_in[k].x[0],
+ p->gparts_in[k].x[1], p->gparts_in[k].x[2], p->nodeID);
+#endif /* SWIFT_DEBUG_CHECKS */
}
if (p->nr_sparts_in > 0) {
if (MPI_Irecv(p->sparts_in, p->nr_sparts_in, spart_mpi_type, p->nodeID,
p->nodeID * proxy_tag_shift + proxy_tag_sparts,
MPI_COMM_WORLD, &p->req_sparts_in) != MPI_SUCCESS)
error("Failed to irecv spart data.");
- // message( "irecv spart data (%i) from node %i." , p->nr_sparts_in ,
- // p->nodeID ); fflush(stdout);
+#ifdef SWIFT_DEBUG_CHECKS
+ message("irecv spart data (%i) from node %i.", p->nr_sparts_in, p->nodeID);
+ fflush(stdout);
+ for (int k = 0; k < p->nr_sparts_in; k++)
+ message(
+ "receiving sparts %lli, x=[%.3e %.3e %.3e], h=%.3e, from node %i.",
+ p->sparts_in[k].id, p->sparts_in[k].x[0], p->sparts_in[k].x[1],
+ p->sparts_in[k].x[2], p->sparts_in[k].h, p->nodeID);
+#endif
}
if (p->nr_bparts_in > 0) {
if (MPI_Irecv(p->bparts_in, p->nr_bparts_in, bpart_mpi_type, p->nodeID,
p->nodeID * proxy_tag_shift + proxy_tag_bparts,
MPI_COMM_WORLD, &p->req_bparts_in) != MPI_SUCCESS)
error("Failed to irecv bpart data.");
- // message( "irecv bpart data (%i) from node %i." , p->nr_bparts_in ,
- // p->nodeID ); fflush(stdout);
+#ifdef SWIFT_DEBUG_CHECKS
+ message("irecv bpart data (%i) from node %i.", p->nr_bparts_in, p->nodeID);
+ fflush(stdout);
+ for (int k = 0; k < p->nr_bparts_in; k++)
+ message(
+ "receiving bparts %lli, x=[%.3e %.3e %.3e], h=%.3e, from node %i.",
+ p->bparts_in[k].id, p->bparts_in[k].x[0], p->bparts_in[k].x[1],
+ p->bparts_in[k].x[2], p->bparts_in[k].h, p->nodeID);
+#endif /* SWIFT_DEBUG_CHECKS */
+ }
+ if (p->nr_sinks_in > 0) {
+ if (MPI_Irecv(p->sinks_in, p->nr_sinks_in, sink_mpi_type, p->nodeID,
+ p->nodeID * proxy_tag_shift + proxy_tag_sinks, MPI_COMM_WORLD,
+ &p->req_sinks_in) != MPI_SUCCESS)
+ error("Failed to irecv sink data.");
+#ifdef SWIFT_DEBUG_CHECKS
+ message("irecv sink data (%i) from node %i.", p->nr_sinks_in, p->nodeID);
+ fflush(stdout);
+ for (int k = 0; k < p->nr_sinks_in; k++)
+ message("receiving sinks %lli, x=[%.3e %.3e %.3e], h=%.3e, from node %i.",
+ p->sinks_in[k].id, p->sinks_in[k].x[0], p->sinks_in[k].x[1],
+ p->sinks_in[k].x[2], p->sinks_in[k].r_cut, p->nodeID);
+#endif /* SWIFT_DEBUG_CHECKS */
}
#else
@@ -987,6 +1091,36 @@ void proxy_bparts_load(struct proxy *p, const struct bpart *bparts, int N) {
p->nr_bparts_out += N;
}
+/**
+ * @brief Load sinks onto a proxy for exchange.
+ *
+ * @param p The #proxy.
+ * @param sinks Pointer to an array of #sink to send.
+ * @param N The number of sinks.
+ */
+void proxy_sinks_load(struct proxy *p, const struct sink *sinks, int N) {
+
+ /* Is there enough space in the buffer? */
+ if (p->nr_sinks_out + N > p->size_sinks_out) {
+ do {
+ p->size_sinks_out *= proxy_buffgrow;
+ } while (p->nr_sinks_out + N > p->size_sinks_out);
+ struct sink *tp;
+ if ((tp = (struct sink *)swift_malloc(
+ "sinks_out", sizeof(struct sink) * p->size_sinks_out)) == NULL)
+ error("Failed to re-allocate sinks_out buffers.");
+ memcpy(tp, p->sinks_out, sizeof(struct sink) * p->nr_sinks_out);
+ swift_free("sinks_out", p->sinks_out);
+ p->sinks_out = tp;
+ }
+
+ /* Copy the parts and xparts data to the buffer. */
+ memcpy(&p->sinks_out[p->nr_sinks_out], sinks, sizeof(struct sink) * N);
+
+ /* Increase the counters. */
+ p->nr_sinks_out += N;
+}
+
/**
* @brief Frees the memory allocated for the particle proxies and sets their
* size back to the initial state.
@@ -1062,6 +1196,21 @@ void proxy_free_particle_buffers(struct proxy *p) {
"bparts_in", sizeof(struct bpart) * p->size_bparts_in)) == NULL)
error("Failed to allocate bparts_in buffers.");
}
+
+ if (p->size_sinks_out > proxy_buffinit) {
+ swift_free("sinks_out", p->sinks_out);
+ p->size_sinks_out = proxy_buffinit;
+ if ((p->sinks_out = (struct sink *)swift_malloc(
+ "sinks_out", sizeof(struct sink) * p->size_sinks_out)) == NULL)
+ error("Failed to allocate sinks_out buffers.");
+ }
+ if (p->size_sinks_in > proxy_buffinit) {
+ swift_free("sinks_in", p->sinks_in);
+ p->size_sinks_in = proxy_buffinit;
+ if ((p->sinks_in = (struct sink *)swift_malloc(
+ "sinks_in", sizeof(struct sink) * p->size_sinks_in)) == NULL)
+ error("Failed to allocate sinks_in buffers.");
+ }
}
/**
@@ -1166,6 +1315,22 @@ void proxy_init(struct proxy *p, int mynodeID, int nodeID) {
error("Failed to allocate bparts_out buffers.");
}
p->nr_bparts_out = 0;
+
+ /* Allocate the sinks send and receive buffers, if needed. */
+ if (p->sinks_in == NULL) {
+ p->size_sinks_in = proxy_buffinit;
+ if ((p->sinks_in = (struct sink *)swift_malloc(
+ "sinks_in", sizeof(struct sink) * p->size_sinks_in)) == NULL)
+ error("Failed to allocate sinks_in buffers.");
+ }
+ p->nr_sinks_in = 0;
+ if (p->sinks_out == NULL) {
+ p->size_sinks_out = proxy_buffinit;
+ if ((p->sinks_out = (struct sink *)swift_malloc(
+ "sinks_out", sizeof(struct sink) * p->size_sinks_out)) == NULL)
+ error("Failed to allocate sinks_out buffers.");
+ }
+ p->nr_sinks_out = 0;
}
/**
@@ -1184,11 +1349,13 @@ void proxy_clean(struct proxy *p) {
swift_free("gparts_out", p->gparts_out);
swift_free("sparts_out", p->sparts_out);
swift_free("bparts_out", p->bparts_out);
+ swift_free("sinks_out", p->sinks_out);
swift_free("parts_in", p->parts_in);
swift_free("xparts_in", p->xparts_in);
swift_free("gparts_in", p->gparts_in);
swift_free("sparts_in", p->sparts_in);
swift_free("bparts_in", p->bparts_in);
+ swift_free("sinks_in", p->sinks_in);
}
/**
diff --git a/src/proxy.h b/src/proxy.h
index d0d2bfb582828eefbbe56d5e7eeba850fa340ffc..46a322e3a49ff53f796353e90aa84852f855acd0 100644
--- a/src/proxy.h
+++ b/src/proxy.h
@@ -28,15 +28,21 @@
#define proxy_buffgrow 1.5
#define proxy_buffinit 100
+/* Number of particle types to exchange with proxies in
+ proxy_parts_exchange_first(). We have parts, gparts, sparts, bparts and
+ sinks to exchange, hence 5 types. */
+#define PROXY_EXCHANGE_NUMBER_PARTICLE_TYPES 5
+
/* Proxy tag arithmetic. */
-#define proxy_tag_shift 8
+#define proxy_tag_shift 9
#define proxy_tag_count 0
#define proxy_tag_parts 1
#define proxy_tag_xparts 2
#define proxy_tag_gparts 3
#define proxy_tag_sparts 4
#define proxy_tag_bparts 5
-#define proxy_tag_cells 6
+#define proxy_tag_sinks 6
+#define proxy_tag_cells 7
/**
* @brief The different reasons a cell can be in a proxy
@@ -71,6 +77,7 @@ struct proxy {
struct gpart *gparts_in, *gparts_out;
struct spart *sparts_in, *sparts_out;
struct bpart *bparts_in, *bparts_out;
+ struct sink *sinks_in, *sinks_out;
int size_parts_in, size_parts_out;
int nr_parts_in, nr_parts_out;
int size_gparts_in, size_gparts_out;
@@ -79,9 +86,12 @@ struct proxy {
int nr_sparts_in, nr_sparts_out;
int size_bparts_in, size_bparts_out;
int nr_bparts_in, nr_bparts_out;
+ int size_sinks_in, size_sinks_out;
+ int nr_sinks_in, nr_sinks_out;
/* Buffer to hold the incomming/outgoing particle counts. */
- int buff_out[4], buff_in[4];
+ int buff_out[PROXY_EXCHANGE_NUMBER_PARTICLE_TYPES],
+ buff_in[PROXY_EXCHANGE_NUMBER_PARTICLE_TYPES];
/* MPI request handles. */
#ifdef WITH_MPI
@@ -91,6 +101,7 @@ struct proxy {
MPI_Request req_gparts_out, req_gparts_in;
MPI_Request req_sparts_out, req_sparts_in;
MPI_Request req_bparts_out, req_bparts_in;
+ MPI_Request req_sinks_out, req_sinks_in;
MPI_Request req_cells_count_out, req_cells_count_in;
MPI_Request req_cells_out, req_cells_in;
#endif
@@ -104,6 +115,7 @@ void proxy_parts_load(struct proxy *p, const struct part *parts,
void proxy_gparts_load(struct proxy *p, const struct gpart *gparts, int N);
void proxy_sparts_load(struct proxy *p, const struct spart *sparts, int N);
void proxy_bparts_load(struct proxy *p, const struct bpart *bparts, int N);
+void proxy_sinks_load(struct proxy *p, const struct sink *sinks, int N);
void proxy_free_particle_buffers(struct proxy *p);
void proxy_parts_exchange_first(struct proxy *p);
void proxy_parts_exchange_second(struct proxy *p);
diff --git a/src/runner_doiact_functions_sinks.h b/src/runner_doiact_functions_sinks.h
index c21c17812849118d1630aeb894591da8c88fdbfc..6ff95ba1deefd310b28f4e24d30bc62d23b5663c 100644
--- a/src/runner_doiact_functions_sinks.h
+++ b/src/runner_doiact_functions_sinks.h
@@ -346,6 +346,10 @@ void DOPAIR1_SINKS_NAIVE(struct runner *r, struct cell *restrict ci,
*/
void DOSELF1_BRANCH_SINKS(struct runner *r, struct cell *c) {
+#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION
+ return;
+#endif
+
const struct engine *restrict e = r->e;
/* Anything to do here? */
@@ -371,6 +375,10 @@ void DOSELF1_BRANCH_SINKS(struct runner *r, struct cell *c) {
*/
void DOPAIR1_BRANCH_SINKS(struct runner *r, struct cell *ci, struct cell *cj) {
+#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION
+ return;
+#endif
+
const struct engine *restrict e = r->e;
const int ci_active = cell_is_active_sinks(ci, e);
@@ -417,6 +425,9 @@ void DOPAIR1_BRANCH_SINKS(struct runner *r, struct cell *ci, struct cell *cj) {
*/
void DOSUB_PAIR1_SINKS(struct runner *r, struct cell *ci, struct cell *cj,
int timer) {
+#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION
+ return;
+#endif
TIMER_TIC;
@@ -496,6 +507,9 @@ void DOSUB_PAIR1_SINKS(struct runner *r, struct cell *ci, struct cell *cj,
* @param gettimer Do we have a timer ?
*/
void DOSUB_SELF1_SINKS(struct runner *r, struct cell *ci, int timer) {
+#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION
+ return;
+#endif
TIMER_TIC;
diff --git a/src/runner_others.c b/src/runner_others.c
index 86b2d6fcbade1dca24eefd5629fb11b333fb8772..82d58f22c46ec5aa0f6e8240e5bd2ef477981da0 100644
--- a/src/runner_others.c
+++ b/src/runner_others.c
@@ -198,6 +198,9 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
*/
void runner_do_star_formation_sink(struct runner *r, struct cell *c,
int timer) {
+#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION
+ return;
+#endif
struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
@@ -567,6 +570,10 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
*/
void runner_do_sink_formation(struct runner *r, struct cell *c) {
+#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION
+ return;
+#endif
+
struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
const struct sink_props *sink_props = e->sink_properties;
diff --git a/src/runner_sinks.c b/src/runner_sinks.c
index fa614f2cbaa23217e2b61f704dbf31fc3039d0f5..d4124c789588bd4b2e56dfe35e0fd24893908f74 100644
--- a/src/runner_sinks.c
+++ b/src/runner_sinks.c
@@ -210,6 +210,9 @@ void runner_do_sinks_gas_swallow(struct runner *r, struct cell *c, int timer) {
*/
void runner_do_sinks_gas_swallow_self(struct runner *r, struct cell *c,
int timer) {
+#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION
+ return;
+#endif
#ifdef SWIFT_DEBUG_CHECKS
if (c->nodeID != r->e->nodeID) error("Running self task on foreign node");
@@ -230,6 +233,9 @@ void runner_do_sinks_gas_swallow_self(struct runner *r, struct cell *c,
*/
void runner_do_sinks_gas_swallow_pair(struct runner *r, struct cell *ci,
struct cell *cj, int timer) {
+#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION
+ return;
+#endif
const struct engine *e = r->e;
@@ -399,6 +405,9 @@ void runner_do_sinks_sink_swallow(struct runner *r, struct cell *c, int timer) {
*/
void runner_do_sinks_sink_swallow_self(struct runner *r, struct cell *c,
int timer) {
+#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION
+ return;
+#endif
#ifdef SWIFT_DEBUG_CHECKS
if (c->nodeID != r->e->nodeID) error("Running self task on foreign node");
@@ -419,6 +428,9 @@ void runner_do_sinks_sink_swallow_self(struct runner *r, struct cell *c,
*/
void runner_do_sinks_sink_swallow_pair(struct runner *r, struct cell *ci,
struct cell *cj, int timer) {
+#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION
+ return;
+#endif
const struct engine *e = r->e;
diff --git a/src/sink/GEAR/sink.h b/src/sink/GEAR/sink.h
index 51ec1b704f75e5510e9b02f6ad6977d40a7f929c..de6bf3636a3a756a5bbf31f35556920ef0c585b4 100644
--- a/src/sink/GEAR/sink.h
+++ b/src/sink/GEAR/sink.h
@@ -172,7 +172,15 @@ __attribute__((always_inline)) INLINE static void sink_init_part(
cpd->E_rot_neighbours[0] = 0.f;
cpd->E_rot_neighbours[1] = 0.f;
cpd->E_rot_neighbours[2] = 0.f;
- cpd->potential = 0.f;
+
+ /* Do not reset the potential to 0. Keep the value computed at the end of the
+ last step. This value is used in runner_iact_nonsym_sink() and
+ runner_iact_sink() to check which particle is at a potential minimum. If you
+ set this value to 0, then we break the check. This value is used instead of
+ gpart->potential because:
+ 1) cpd->potential does not break MPI, while gpart->potential does
+ 2) gpart->potential is not yet computed in runner_iact_X_sink(). */
+ /* cpd->potential = 0.f; */
cpd->E_mec_bound = 0.f; /* Gravitationally bound particles will have
E_mec_bound < 0. This is checked before comparing
any other value with this one. So no need to put
diff --git a/src/sink/GEAR/sink_iact.h b/src/sink/GEAR/sink_iact.h
index fc75ba8842bb7ceaae6d789b277a76b21206c29f..fe3791cf964155d0ae7b34276cec09dc01f5a967 100644
--- a/src/sink/GEAR/sink_iact.h
+++ b/src/sink/GEAR/sink_iact.h
@@ -33,8 +33,6 @@
* In GEAR: This function deactivates the sink formation ability of #part not
* at a potential minimum.
*
- * Note: This functions breaks MPI.
- *
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
@@ -58,12 +56,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_sink(
/* if the distance is less than the cut off radius */
if (r < cut_off_radius) {
-
- /*
- * NOTE: Those lines break MPI
- */
- float potential_i = pi->gpart->potential;
- float potential_j = pj->gpart->potential;
+ float potential_i = pi->sink_data.potential;
+ float potential_j = pj->sink_data.potential;
/* prevent the particle with the largest potential to form a sink */
if (potential_i > potential_j) {
@@ -107,9 +101,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink(
const float r = sqrtf(r2);
if (r < cut_off_radius) {
-
- float potential_i = pi->gpart->potential;
- float potential_j = pj->gpart->potential;
+ float potential_i = pi->sink_data.potential;
+ float potential_j = pj->sink_data.potential;
/* if the potential is larger
* prevent the particle to form a sink */
@@ -155,6 +148,9 @@ runner_iact_nonsym_sinks_gas_density(
*
* Note: Energies are computed with physical quantities, not the comoving ones.
*
+ * MPI note: This functions invokes the gpart. Hence, it must be performed only
+ * on the local node (similarly to runner_iact_nonsym_bh_bh_repos()).
+ *
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param ri Comoving cut off radius of particle i.
@@ -286,6 +282,9 @@ runner_iact_nonsym_sinks_sink_swallow(
*
* Note: Energies are computed with physical quantities, not the comoving ones.
*
+ * MPI note: This functions invokes the gpart. Hence, it must be performed only
+ * on the local node (similarly to runner_iact_nonsym_bh_gas_repos()).
+ *
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param ri Comoving cut off radius of particle i.
diff --git a/src/space.c b/src/space.c
index 72125200b25aa0e01edd66aeed508c89e325caf9..95d029c73da5828759b39a971fb1dee12cbef0c2 100644
--- a/src/space.c
+++ b/src/space.c
@@ -150,6 +150,11 @@ void space_free_foreign_parts(struct space *s, const int clear_cell_pointers) {
s->size_bparts_foreign = 0;
s->bparts_foreign = NULL;
}
+ if (s->sinks_foreign != NULL) {
+ swift_free("sinks_foreign", s->sinks_foreign);
+ s->size_sinks_foreign = 0;
+ s->sinks_foreign = NULL;
+ }
if (clear_cell_pointers) {
for (int k = 0; k < s->e->nr_proxies; k++) {
for (int j = 0; j < s->e->proxies[k].nr_cells_in; j++) {
@@ -2175,6 +2180,7 @@ void space_check_drift_point(struct space *s, integertime_t ti_drift,
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);
+ space_map_cells_pre(s, 1, cell_check_sink_drift_point, &ti_drift);
if (multipole)
space_map_cells_pre(s, 1, cell_check_multipole_drift_point, &ti_drift);
#else
@@ -2321,6 +2327,55 @@ void space_check_bpart_swallow_mapper(void *map_data, int nr_bparts,
#endif
}
+/**
+ * @brief #threadpool mapper function for the swallow debugging check
+ */
+void space_check_part_sink_swallow_mapper(void *map_data, int nr_parts,
+ void *extra_data) {
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Unpack the data */
+ struct part *restrict parts = (struct part *)map_data;
+
+ /* Verify that all particles have been swallowed or are untouched */
+ for (int k = 0; k < nr_parts; k++) {
+
+ if (parts[k].time_bin == time_bin_inhibited) continue;
+
+ const long long swallow_id = sink_get_part_swallow_id(&parts[k].sink_data);
+
+ if (swallow_id != -1)
+ error("Particle has not been swallowed! id=%lld", parts[k].id);
+ }
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
+/**
+ * @brief #threadpool mapper function for the swallow debugging check
+ */
+void space_check_sink_sink_swallow_mapper(void *map_data, int nr_sinks,
+ void *extra_data) {
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Unpack the data */
+ struct sink *restrict sinks = (struct sink *)map_data;
+
+ /* Verify that all particles have been swallowed or are untouched */
+ for (int k = 0; k < nr_sinks; k++) {
+
+ if (sinks[k].time_bin == time_bin_inhibited) continue;
+
+ const long long swallow_id =
+ sink_get_sink_swallow_id(&sinks[k].merger_data);
+
+ if (swallow_id != -1)
+ error("Sink particle has not been swallowed! id=%lld", sinks[k].id);
+ }
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
/**
* @brief Checks that all particles have their swallow flag in a "no swallow"
* state.
@@ -2339,6 +2394,16 @@ void space_check_swallow(struct space *s) {
threadpool_map(&s->e->threadpool, space_check_bpart_swallow_mapper, s->bparts,
s->nr_bparts, sizeof(struct bpart), threadpool_auto_chunk_size,
/*extra_data=*/NULL);
+
+ threadpool_map(&s->e->threadpool, space_check_part_sink_swallow_mapper,
+ s->parts, s->nr_parts, sizeof(struct part),
+ threadpool_auto_chunk_size,
+ /*extra_data=*/NULL);
+
+ threadpool_map(&s->e->threadpool, space_check_sink_sink_swallow_mapper,
+ s->sinks, s->nr_sinks, sizeof(struct sink),
+ threadpool_auto_chunk_size,
+ /*extra_data=*/NULL);
#else
error("Calling debugging code without debugging flag activated.");
#endif
@@ -2437,6 +2502,7 @@ void space_clean(struct space *s) {
swift_free("sparts_foreign", s->sparts_foreign);
swift_free("gparts_foreign", s->gparts_foreign);
swift_free("bparts_foreign", s->bparts_foreign);
+ swift_free("sinks_foreign", s->sinks_foreign);
#endif
free(s->cells_sub);
free(s->multipoles_sub);
@@ -2611,6 +2677,8 @@ void space_struct_restore(struct space *s, FILE *stream) {
s->size_sparts_foreign = 0;
s->bparts_foreign = NULL;
s->size_bparts_foreign = 0;
+ s->sinks_foreign = NULL;
+ s->size_sinks_foreign = 0;
#endif
/* More things to read. */
@@ -2727,9 +2795,10 @@ void space_write_cell(const struct space *s, FILE *f, const struct cell *c) {
/* Write line for current cell */
fprintf(f, "%lld,%lld,%i,", c->cellID, parent, c->nodeID);
- fprintf(f, "%i,%i,%i,%s,%s,%g,%g,%g,%g,%g,%g, ", c->hydro.count,
- c->stars.count, c->grav.count, superID, hydro_superID, c->loc[0],
- c->loc[1], c->loc[2], c->width[0], c->width[1], c->width[2]);
+ fprintf(f, "%i,%i,%i,%i,%s,%s,%g,%g,%g,%g,%g,%g, ", c->hydro.count,
+ c->stars.count, c->grav.count, c->sinks.count, superID, hydro_superID,
+ c->loc[0], c->loc[1], c->loc[2], c->width[0], c->width[1],
+ c->width[2]);
fprintf(f, "%g, %g, %i, %i\n", c->hydro.h_max, c->stars.h_max, c->depth,
c->maxdepth);
@@ -2761,7 +2830,7 @@ void space_write_cell_hierarchy(const struct space *s, int j) {
if (engine_rank == 0) {
fprintf(f, "name,parent,mpi_rank,");
fprintf(f,
- "hydro_count,stars_count,gpart_count,super,hydro_super,"
+ "hydro_count,stars_count,gpart_count,sinks_count,super,hydro_super,"
"loc1,loc2,loc3,width1,width2,width3,");
fprintf(f, "hydro_h_max,stars_h_max,depth,maxdepth\n");
diff --git a/src/space.h b/src/space.h
index df58eeda4b4bca6d9909786247877cf1caafdec2..8a85277f82d6d9169f6f984d02671f870b9a2961 100644
--- a/src/space.h
+++ b/src/space.h
@@ -349,6 +349,10 @@ struct space {
struct bpart *bparts_foreign;
size_t nr_bparts_foreign, size_bparts_foreign;
+ /*! Buffers for sink-parts that we will receive from foreign cells. */
+ struct sink *sinks_foreign;
+ size_t nr_sinks_foreign, size_sinks_foreign;
+
#endif
};
diff --git a/src/space_cell_index.c b/src/space_cell_index.c
index c9d7cdd07fdcce9a2333d40e9d94fa10e064a75d..c0ad4f2df00160cc341d585cf6295c2728e4795a 100644
--- a/src/space_cell_index.c
+++ b/src/space_cell_index.c
@@ -692,8 +692,8 @@ void space_sinks_get_cell_index_mapper(void *map_data, int nr_sinks,
if (count_extra_sink) atomic_add(&data->count_extra_sink, count_extra_sink);
/* Write back the minimal part mass and velocity sum */
- atomic_min_f(&s->min_spart_mass, min_mass);
- atomic_add_f(&s->sum_spart_vel_norm, sum_vel_norm);
+ atomic_min_f(&s->min_sink_mass, min_mass);
+ atomic_add_f(&s->sum_sink_vel_norm, sum_vel_norm);
}
/**
diff --git a/src/space_extras.c b/src/space_extras.c
index 9d8a88952e5a3d2a670ba63f96001dc27eda26c7..64ea69647ffc1f88ce600703a931a53626dea3d5 100644
--- a/src/space_extras.c
+++ b/src/space_extras.c
@@ -359,6 +359,7 @@ void space_allocate_extras(struct space *s, int verbose) {
bzero(&s->sinks[i], sizeof(struct sink));
s->sinks[i].time_bin = time_bin_not_created;
s->sinks[i].id = -42;
+ sink_mark_sink_as_not_swallowed(&s->sinks[i].merger_data);
}
/* Put the spare particles in their correct cell */
diff --git a/src/space_first_init.c b/src/space_first_init.c
index fd36583715849c662ac56ae47d28ceff4396074a..3232e445f0fb21fb1dfeae14d91786efa2ed8688 100644
--- a/src/space_first_init.c
+++ b/src/space_first_init.c
@@ -484,6 +484,9 @@ void space_first_init_sinks_mapper(void *restrict map_data, int count,
/* Initialize the sinks */
sink_first_init_sink(&sink[k], props, e);
+ /* And the sink merger markers */
+ sink_mark_sink_as_not_swallowed(&sink[k].merger_data);
+
/* Also initialize the chemistry */
chemistry_first_init_sink(chemistry, &sink[k]);
diff --git a/src/space_rebuild.c b/src/space_rebuild.c
index 36d784c8b34c393d4260fbd54ae5fa02b5c4c771..7d458f6858ea64ea9970a376d177c2211cc5ea0a 100644
--- a/src/space_rebuild.c
+++ b/src/space_rebuild.c
@@ -491,17 +491,20 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
size_t nr_gparts_exchanged = s->nr_gparts - nr_gparts;
size_t nr_sparts_exchanged = s->nr_sparts - nr_sparts;
size_t nr_bparts_exchanged = s->nr_bparts - nr_bparts;
+ size_t nr_sinks_exchanged = s->nr_sinks - nr_sinks;
engine_exchange_strays(s->e, nr_parts, &h_index[nr_parts],
&nr_parts_exchanged, nr_gparts, &g_index[nr_gparts],
&nr_gparts_exchanged, nr_sparts, &s_index[nr_sparts],
&nr_sparts_exchanged, nr_bparts, &b_index[nr_bparts],
- &nr_bparts_exchanged);
+ &nr_bparts_exchanged, nr_sinks,
+ &sink_index[nr_sinks], &nr_sinks_exchanged);
/* Set the new particle counts. */
s->nr_parts = nr_parts + nr_parts_exchanged;
s->nr_gparts = nr_gparts + nr_gparts_exchanged;
s->nr_sparts = nr_sparts + nr_sparts_exchanged;
s->nr_bparts = nr_bparts + nr_bparts_exchanged;
+ s->nr_sinks = nr_sinks + nr_sinks_exchanged;
} else {
#ifdef SWIFT_DEBUG_CHECKS
@@ -511,6 +514,8 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
error("Number of sparts changing after repartition");
if (s->nr_gparts != nr_gparts)
error("Number of gparts changing after repartition");
+ if (s->nr_sinks != nr_sinks)
+ error("Number of sinks changing after repartition");
#endif
}
@@ -521,6 +526,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
cell_spart_counts[k] = 0;
cell_gpart_counts[k] = 0;
cell_bpart_counts[k] = 0;
+ cell_sink_counts[k] = 0;
}
}
@@ -547,16 +553,27 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
/* Re-allocate the index array for the bparts if needed.. */
- if (s->nr_bparts + 1 > s_index_size) {
+ if (s->nr_bparts + 1 > b_index_size) {
int *bind_new;
if ((bind_new = (int *)swift_malloc(
"b_index", sizeof(int) * (s->nr_bparts + 1))) == NULL)
- error("Failed to allocate temporary s-particle indices.");
+ error("Failed to allocate temporary b-particle indices.");
memcpy(bind_new, b_index, sizeof(int) * nr_bparts);
swift_free("b_index", b_index);
b_index = bind_new;
}
+ /* Re-allocate the index array for the sinks if needed.. */
+ if (s->nr_sinks + 1 > sink_index_size) {
+ int *sink_ind_new;
+ if ((sink_ind_new = (int *)swift_malloc(
+ "sink_index", sizeof(int) * (s->nr_sinks + 1))) == NULL)
+ error("Failed to allocate temporary sink-particle indices.");
+ memcpy(sink_ind_new, sink_index, sizeof(int) * nr_sinks);
+ swift_free("sink_index", sink_index);
+ sink_index = sink_ind_new;
+ }
+
const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
const double ih[3] = {s->iwidth[0], s->iwidth[1], s->iwidth[2]};
@@ -602,6 +619,20 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
nr_bparts = s->nr_bparts;
+ /* Assign each received sink to its cell. */
+ for (size_t k = nr_sinks; k < s->nr_sinks; k++) {
+ const struct sink *const sink = &s->sinks[k];
+ sink_index[k] = cell_getid(cdim, sink->x[0] * ih[0], sink->x[1] * ih[1],
+ sink->x[2] * ih[2]);
+ cell_sink_counts[sink_index[k]]++;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (cells_top[sink_index[k]].nodeID != local_nodeID)
+ error("Received sink-part that does not belong to me (nodeID=%i).",
+ cells_top[sink_index[k]].nodeID);
+#endif
+ }
+ nr_sinks = s->nr_sinks;
+
#else /* WITH_MPI */
/* Update the part, spart and bpart counters */
@@ -716,14 +747,14 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
error("Inhibited particle sorted into a cell!");
/* New cell index */
- const int new_bind =
+ const int new_sink_ind =
cell_getid(s->cdim, sink->x[0] * s->iwidth[0],
sink->x[1] * s->iwidth[1], sink->x[2] * s->iwidth[2]);
/* New cell of this sink */
- const struct cell *c = &s->cells_top[new_bind];
+ const struct cell *c = &s->cells_top[new_sink_ind];
- if (sink_index[k] != new_bind)
+ if (sink_index[k] != new_sink_ind)
error("sink's new cell index not matching sorted index.");
if (sink->x[0] < c->loc[0] || sink->x[0] > c->loc[0] + c->width[0] ||
@@ -931,6 +962,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* Store the state at rebuild time */
c->stars.parts_rebuild = c->stars.parts;
+ c->sinks.parts_rebuild = c->sinks.parts;
c->grav.parts_rebuild = c->grav.parts;
c->hydro.count_total = c->hydro.count + space_extra_parts;
diff --git a/src/space_recycle.c b/src/space_recycle.c
index 1a7a42c830e3636815feea2b4ad5e1d70bdccb45..b510eb9dced0f3ee4e9eee47a0d8c434962e0c03 100644
--- a/src/space_recycle.c
+++ b/src/space_recycle.c
@@ -193,6 +193,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->grav.parts = NULL;
c->grav.parts_rebuild = NULL;
c->sinks.parts = NULL;
+ c->sinks.parts_rebuild = NULL;
c->stars.parts = NULL;
c->stars.parts_rebuild = NULL;
c->black_holes.parts = NULL;
diff --git a/src/space_regrid.c b/src/space_regrid.c
index 487fe7c0e38495ac85622156615308cd4179ec01..dc86468e335b3cc2e50d6250327f532bc64e22eb 100644
--- a/src/space_regrid.c
+++ b/src/space_regrid.c
@@ -64,7 +64,8 @@ void space_regrid(struct space *s, int verbose) {
if (c->black_holes.h_max > h_max) {
h_max = c->black_holes.h_max;
}
- if (c->sinks.r_cut_max > h_max) {
+ /* Note: For sinks, r_cut <=> h*kernel_gamma */
+ if (c->sinks.r_cut_max > h_max * kernel_gamma) {
h_max = c->sinks.r_cut_max / kernel_gamma;
}
}
@@ -82,7 +83,9 @@ void space_regrid(struct space *s, int verbose) {
if (c->nodeID == engine_rank && c->black_holes.h_max > h_max) {
h_max = c->black_holes.h_max;
}
- if (c->nodeID == engine_rank && c->sinks.r_cut_max > h_max) {
+ /* Note: For sinks, r_cut <=> h*kernel_gamma */
+ if (c->nodeID == engine_rank &&
+ c->sinks.r_cut_max > h_max * kernel_gamma) {
h_max = c->sinks.r_cut_max / kernel_gamma;
}
}
diff --git a/src/space_unique_id.c b/src/space_unique_id.c
index bfcb733cc6daa698e2744251d7f76b5dc93068e0..481d8a362799e310422e6020361d8ea7fc1ab913 100644
--- a/src/space_unique_id.c
+++ b/src/space_unique_id.c
@@ -82,9 +82,10 @@ void space_update_unique_id(struct space *s) {
#endif // WITH_MPI
/* Compute the size of each batch. */
- const long long batch_size = (space_extra_parts + space_extra_sparts +
- space_extra_gparts + space_extra_bparts) *
- s->nr_cells;
+ const long long batch_size =
+ (space_extra_parts + space_extra_sparts + space_extra_gparts +
+ space_extra_bparts + space_extra_sinks) *
+ s->nr_cells;
/* Get a new batch. */
if (require_new_batch) {
@@ -210,9 +211,10 @@ void space_init_unique_id(struct space *s, int nr_nodes) {
s->unique_id.global_next_id++;
/* Compute the size of each batch. */
- const long long batch_size = (space_extra_parts + space_extra_sparts +
- space_extra_gparts + space_extra_bparts) *
- s->nr_cells;
+ const long long batch_size =
+ (space_extra_parts + space_extra_sparts + space_extra_gparts +
+ space_extra_bparts + space_extra_sinks) *
+ s->nr_cells;
/* Compute the initial batchs (each rank has 2 batchs). */
if (s->unique_id.global_next_id > LLONG_MAX - 2 * engine_rank * batch_size) {
diff --git a/swift.c b/swift.c
index b63941cd63538e3226669bfed06bdb1837c73411..41c96e0afbc4ca5ad2a119e2089397a5fcbd540e 100644
--- a/swift.c
+++ b/swift.c
@@ -528,11 +528,17 @@ int main(int argc, char *argv[]) {
#endif
#ifdef WITH_MPI
+#ifdef SWIFT_DEBUG_CHECKS
+ if (with_sinks) {
+ pretime_message("Warning: sink particles are are WIP yet with MPI.");
+ }
+#else
if (with_sinks) {
pretime_message("Error: sink particles are not available yet with MPI.");
return 1;
}
-#endif
+#endif /* SWIFT_DEBUG_CHECKS */
+#endif /* WITH_MPI */
if (with_sinks && with_star_formation) {
pretime_message(
diff --git a/tests/difffloat.py b/tests/difffloat.py
index cfa16970730bc365b284a78f8a757bd459a9b741..6f60a66535b39e9aa4d06951cc088f975435d46a 100644
--- a/tests/difffloat.py
+++ b/tests/difffloat.py
@@ -18,6 +18,7 @@
##############################################################################
from numpy import *
+
del min
import sys
diff --git a/theory/SinkParticles/plot_imf_image.py b/theory/SinkParticles/plot_imf_image.py
index 25da38cab2b804a651c5f745063e6ba101b38a52..bc1da9cccf1c3fe60af2661e4cf65a96fb949338 100644
--- a/theory/SinkParticles/plot_imf_image.py
+++ b/theory/SinkParticles/plot_imf_image.py
@@ -14,18 +14,17 @@ import matplotlib.pyplot as plt
import matplotlib
# For swift doc, use the following and not the styleheet
-plt.rcParams['text.usetex'] = True
-plt.rcParams['font.family'] = 'serif'
-plt.rcParams['font.serif'] = ['Computer Modern Roman']
+plt.rcParams["text.usetex"] = True
+plt.rcParams["font.family"] = "serif"
+plt.rcParams["font.serif"] = ["Computer Modern Roman"]
mmin = 0.5
mmax = 300
-matplotlib.rcParams.update({'font.size': 16})
+matplotlib.rcParams.update({"font.size": 16})
figsize = (6.4, 4.8)
-fig, ax = plt.subplots(num=1, ncols=1, nrows=1,
- figsize=figsize, layout="tight")
+fig, ax = plt.subplots(num=1, ncols=1, nrows=1, figsize=figsize, layout="tight")
ax.set_xlim([0.3, 400])
ax.set_ylim([1, 5e4])
@@ -34,34 +33,34 @@ ax.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax.set_xlabel("$M_{\star}$ $[M_\odot]$")
ax.set_ylabel("d$N/$d$M$ [arbitrary units]")
-ax.set_xscale('log')
-ax.set_yscale('log')
+ax.set_xscale("log")
+ax.set_yscale("log")
# theoretical imf
s = -1.3
-bins = 10**np.linspace(np.log10(mmin), np.log10(mmax), 100)
-n = 0.9*10000*bins**s
+bins = 10 ** np.linspace(np.log10(mmin), np.log10(mmax), 100)
+n = 0.9 * 10000 * bins ** s
ax.plot(bins, n, "k--")
-bins = 10**np.linspace(np.log10(mmin), np.log10(8), 100)
-n = 0.9*10000*bins**s
+bins = 10 ** np.linspace(np.log10(mmin), np.log10(8), 100)
+n = 0.9 * 10000 * bins ** s
ax.fill_between(bins, 0.1, n, color="red", alpha=0.1)
-bins = 10**np.linspace(np.log10(8), np.log10(mmax), 100)
-n = 0.9*10000*bins**s
+bins = 10 ** np.linspace(np.log10(8), np.log10(mmax), 100)
+n = 0.9 * 10000 * bins ** s
ax.fill_between(bins, 0.1, n, color="blue", alpha=0.1)
-ax.text(2, 1e2, r"$M_{\rm c}$", horizontalalignment='center')
-ax.text(50, 2, r"$M_{\rm d}$", horizontalalignment='center')
+ax.text(2, 1e2, r"$M_{\rm c}$", horizontalalignment="center")
+ax.text(50, 2, r"$M_{\rm d}$", horizontalalignment="center")
# Add limit
-ax.vlines(x=8, ymin=0, ymax=600, color='k', linestyle='-')
+ax.vlines(x=8, ymin=0, ymax=600, color="k", linestyle="-")
# Add text to the vertical line
-ax.text(8, 800, r"$m_{t}$", horizontalalignment='center')
+ax.text(8, 800, r"$m_{t}$", horizontalalignment="center")
# fig.patch.set_facecolor('none') # Remove figure background
# ax.set_facecolor('none') # Remove axes background
-plt.savefig('sink_imf.png', dpi=300, bbox_inches='tight')
+plt.savefig("sink_imf.png", dpi=300, bbox_inches="tight")
plt.close()
diff --git a/tools/data/cell_hierarchy.html b/tools/data/cell_hierarchy.html
index 6c9bcdb995be7970948cf9b9544ae0a59d41caf0..051451091d70d7e6bcb000b44b71f73c2779caa6 100644
--- a/tools/data/cell_hierarchy.html
+++ b/tools/data/cell_hierarchy.html
@@ -162,6 +162,7 @@
"MPI rank:" + n.mpi_rank + "<br/>" +
"Part: " + n.hydro_count + "<br/>" +
"Spart: " + n.stars_count + "<br/>" +
+ "Sink: " + n.sinks_count + "<br/>" +
"Super: " + n.super + "<br/>" +
"Super Hydro: " + n.hydro_super + "<br/>" +
"Loc: " + n.loc1 + ", " + n.loc2 + ", " + n.loc3 + "<br/>" +
diff --git a/tools/plot_task_dependencies.py b/tools/plot_task_dependencies.py
index 17013ca2f131ddf648d81f4fa70158df65bbd610..fa9fdb7fd1e5a04eb340f9fee0801918a6f3370c 100755
--- a/tools/plot_task_dependencies.py
+++ b/tools/plot_task_dependencies.py
@@ -208,7 +208,12 @@ def task_is_hydro(name):
return True
if "_part" in name:
return True
- if "density" in name and "stars" not in name and "sink" not in name and "bh" not in name:
+ if (
+ "density" in name
+ and "stars" not in name
+ and "sink" not in name
+ and "bh" not in name
+ ):
return True
if "rho" in name and "bpart" not in name:
return True