diff --git a/doc/RTD/source/FriendsOfFriends/algorithm_description.rst b/doc/RTD/source/FriendsOfFriends/algorithm_description.rst
index bf7fd77fbc36f026088c317bf9b50c24c1406061..17a667b48be48a4e4dd6d0be5c71fbf306933a3d 100644
--- a/doc/RTD/source/FriendsOfFriends/algorithm_description.rst
+++ b/doc/RTD/source/FriendsOfFriends/algorithm_description.rst
@@ -19,8 +19,19 @@ friends (its *friends-of-friends*). This creates networks of linked particles
which are called *groups*. The size (or length) of
a group is the number of particles in that group. If a particle does not
find any other particle within ``l`` then it forms its own group of
-size 1. For a given distribution of particles the resulting list of
-groups is unique and unambiguously defined.
+size 1. **For a given distribution of particles the resulting list of
+groups is unique and unambiguously defined.**
+
+In our implementation, we use three separate categories influencing their
+behaviour in the FOF code:
+ * ``linkable`` particles which behave as described above.
+ * ``attachable`` particles which can `only` form a link with the `nearest` ``linkable`` particle they find.
+ * And the others which are ignored entirely.
+
+The category of each particle type is specified at run time in the parameter
+file. The classic scenario for the two categories is to run FOF on the dark
+matter particles (i.e. they are `linkable`) and then attach the gas, stars and
+black holes to their nearest DM (i.e. the baryons are `attachable`).
Small groups are typically discarded, the final catalogue only contains
objects with a length above a minimal threshold, typically of the
@@ -36,20 +47,25 @@ domain decomposition and tree structure that is created for the other
parts of the code. The tree can be easily used to find neighbours of
particles within the linking length.
-Depending on the application, the choice of linking length and
-minimal group size can vary. For cosmological applications, bound
-structures (dark matter haloes) are traditionally identified using a
-linking length expressed as :math:`0.2` of the mean inter-particle
-separation :math:`d` in the simulation which is given by :math:`d =
-\sqrt[3]{\frac{V}{N}}`, where :math:`N` is the number of particles in
-the simulation and :math:`V` is the simulation (co-moving)
-volume. Usually only dark matter particles are considered for the
-number :math:`N`. Other particle types are linked but do not
-participate in the calculation of the linking length. Experience shows
-that this produces groups that are similar to the commonly adopted
-(but much more complex) definition of virialised haloes. A minimal
-group length of :math:`32` is often adopted in order to get a robust
-catalogue of haloes and compute a good halo mass function.
+Depending on the application, the choice of linking length and minimal group
+size can vary. For cosmological applications, bound structures (dark matter
+haloes) are traditionally identified using a linking length expressed as
+:math:`0.2` of the mean inter-particle separation :math:`d` in the simulation
+which is given by :math:`d = \sqrt[3]{\frac{V}{N}}`, where :math:`N` is the
+number of particles in the simulation and :math:`V` is the simulation
+(co-moving) volume. Experience shows that this produces groups that are similar
+to the commonly adopted (but much more complex) definition of virialised
+haloes. A minimal group length of :math:`32` is often adopted in order to get a
+robust catalogue of haloes and compute a good halo mass function. Usually only
+dark matter particles are considered for the number :math:`N`. In practice, the
+mean inter-particle separation is evaluated based on the cosmology adopted in
+the simulation. We use: :math:`d=\sqrt[3]{\frac{m_{\rm DM}}{\Omega_{\rm cdm}
+\rho_{\rm crit}}}` for simulations with baryonic particles and
+:math:`d=\sqrt[3]{\frac{m_{\rm DM}}{(\Omega_{\rm cdm} + \Omega_{\rm b})
+\rho_{\rm crit}}}` for DMO simulations. In both cases, :math:`m_{\rm DM}` is the
+mean mass of the DM particles. Using this definition (rather than basing in on
+:math:`N`) makes the code robust to zoom-in scenarios where the entire volume is
+not filled with particles.
For non-cosmological applications of the FOF algorithm, the choice of
the linking length is more difficult and left to the user. The choice
diff --git a/doc/RTD/source/FriendsOfFriends/on_the_fly_fof.rst b/doc/RTD/source/FriendsOfFriends/on_the_fly_fof.rst
index 9961cca88366e398bb05dace29633a7293a6ee77..06eca6769bcbaf993dcbdd50e565e2c25d79d4c7 100644
--- a/doc/RTD/source/FriendsOfFriends/on_the_fly_fof.rst
+++ b/doc/RTD/source/FriendsOfFriends/on_the_fly_fof.rst
@@ -10,8 +10,9 @@ The main purpose of the on-the-fly FOF is to identify haloes during a
cosmological simulation in order to seed some of them with black holes
based on physical considerations.
-**In this mode, no group catalogue is written to the disk. The resulting list
-of haloes is only used internally by SWIFT.**
+.. warning::
+ In this mode, no group catalogue is written to the disk. The resulting list
+ of haloes is only used internally by SWIFT.
Note that a catalogue can nevertheless be written after every seeding call by
setting the optional parameter ``dump_catalogue_when_seeding``.
diff --git a/doc/RTD/source/FriendsOfFriends/parameter_description.rst b/doc/RTD/source/FriendsOfFriends/parameter_description.rst
index d9820a74d2bea768904b5fdeb35aacd01e04efd1..0e2afa97ceb7806d2851fe68ed4fc8e9d22ae082 100644
--- a/doc/RTD/source/FriendsOfFriends/parameter_description.rst
+++ b/doc/RTD/source/FriendsOfFriends/parameter_description.rst
@@ -20,6 +20,12 @@ absolute value using the parameter ``absolute_linking_length``. This is
expressed in internal units. This value will be ignored (and the ratio of
the mean inter-particle separation will be used) when set to ``-1``.
+The categories of particles are specified using the ``linking_types`` and
+``attaching_types`` arrays. They are of the length of the number of particle
+types in SWIFT (currently 7) and specify for each type using ``1`` or ``0``
+whether or not the given particle type is in this category. Types not present
+in either category are ignored entirely.
+
The second important parameter is the minimal size of groups to retain in
the catalogues. This is given in terms of number of particles (of all types)
via the parameter ``min_group_size``. When analysing simulations, to
@@ -98,10 +104,12 @@ A full FOF section of the YAML parameter file looks like:
time_first: 0.2 # Time of first FoF black hole seeding calls.
delta_time: 1.005 # Time between consecutive FoF black hole seeding calls.
min_group_size: 256 # The minimum no. of particles required for a group.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Which particle types to consider for linking (here only DM)
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Which particle types to consider for attaching (here gas, stars, and BHs)
linking_length_ratio: 0.2 # Linking length in units of the main inter-particle separation.
+ seed_black_holes_enabled: 0 # Do not seed black holes when running FOF
black_hole_seed_halo_mass_Msun: 1.5e10 # Minimal halo mass in which to seed a black hole (in solar masses).
dump_catalogue_when_seeding: 0 # (Optional) Write a FOF catalogue when seeding black holes. Defaults to 0 if unspecified.
absolute_linking_length: -1. # (Optional) Absolute linking length (in internal units).
group_id_default: 2147483647 # (Optional) Sets the group ID of particles in groups below the minimum size.
group_id_offset: 1 # (Optional) Sets the offset of group ID labelling. Defaults to 1 if unspecified.
- seed_black_holes_enabled: 0 # Do not seed black holes when running FOF
diff --git a/doc/RTD/source/FriendsOfFriends/stand_alone_fof.rst b/doc/RTD/source/FriendsOfFriends/stand_alone_fof.rst
index d5b341409b1553dcfb4062765d83f7ba14e18570..ecf82634aed909a77c81cf27ea438aaec9b98338 100644
--- a/doc/RTD/source/FriendsOfFriends/stand_alone_fof.rst
+++ b/doc/RTD/source/FriendsOfFriends/stand_alone_fof.rst
@@ -11,17 +11,17 @@ compiled by configuring the code with the option
``--enable-stand-alone-fof``. The ``fof`` and ``fof_mpi`` executables
will then be generated alongside the regular SWIFT ones.
-The executable takes a parameter file as an argument. It will then
-read the snapshot specified in the parameter file and extract all
-the dark matter particles by default. FOF is then run on these
-particles and a catalogue of groups is written to disk. Additional
-particle types can be read and processed by the stand-alone FOF
-code by adding any of the following runtime parameters to the
-command line:
+The executable takes a parameter file as an argument. It will then read the
+snapshot specified in the parameter file (specified as an initial condition
+file) and extract all the dark matter particles by default. FOF is then run on
+these particles and a catalogue of groups is written to disk. Additional
+particle types can be read and processed by the stand-alone FOF code by adding
+any of the following runtime parameters to the command line:
* ``--hydro``: Read and process the gas particles,
* ``--stars``: Read and process the star particles,
* ``--black-holes``: Read and process the black hole particles,
+ * ``--sinks``: Read and process the sink particles,
* ``--cosmology``: Consider cosmological terms.
Running with cosmology is necessary when using a linking length based
@@ -34,3 +34,10 @@ internal units). The FOF code will also write a snapshot with an
additional field for each particle. This contains the ``GroupID`` of
each particle and can be used to find all the particles in a given
halo and to link them to the information stored in the catalogue.
+
+.. warning::
+
+ Note that since not all particle properties are read in stand-alone
+ mode, not all particle properties will be written to the snapshot generated
+ by the stand-alone FOF.
+
diff --git a/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml b/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml
index 0fd2c5f589557b1738631ca43c7f8c46ff91bdb7..85ae233751d7df7e915e0df8fc3cbf75eb3dbc88 100644
--- a/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml
+++ b/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml
@@ -76,6 +76,8 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.0e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.05 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.00751 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
Scheduler:
max_top_level_cells: 64
diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index 2e33e031e1c56423f55760d32c5a39e1444de312..0742866a4ab5be0367a6adbc99aed25ebccf1fe3 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -76,6 +76,8 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.0e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.05 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.00751 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
Scheduler:
max_top_level_cells: 8
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index 77355071563c56be04be67fcb617613ef6cc44af..e45a3dd456d7309d804d51e263d5cf5b79b832b4 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -76,6 +76,8 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.0e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.05 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.00751 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
Scheduler:
max_top_level_cells: 16
diff --git a/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml
index 743fa0d5dede782ecb1445d861d60f6190171037..cdb4d98e80c243ff6b40407ca20eace355d0ae7f 100644
--- a/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml
@@ -76,6 +76,8 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.0e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.05 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.00751 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
Scheduler:
max_top_level_cells: 16
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index 76f454627764c12224c1724892a5d5b05cf65d32..29ca9fe3909d9ff3479af7c876b7bdd65ca678ac 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -76,6 +76,8 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.0e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.05 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.00751 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
Scheduler:
max_top_level_cells: 32
diff --git a/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml
index 97af6a1f60019d1ebd5fa6e2670b36ade68ef572..a14bb9e5ac51729afa55f3c0217c1ebdd1c6fedf 100644
--- a/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml
@@ -76,6 +76,8 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.0e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.05 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.00751 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
Scheduler:
max_top_level_cells: 16
diff --git a/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml b/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml
index 4b8446c30d5401f31c9bc3650b9389fb47fc7c14..1a22107d119eb3389d8898d82c28c6862966a1bb 100644
--- a/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml
@@ -76,6 +76,8 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.0e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.05 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.00751 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
Scheduler:
max_top_level_cells: 8
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index 96fe8efe01082f79a676110c5ae6e59e25f06dc1..032b8b6e7d0c1a5ebc07c71c905ac7bb0787ddac 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -83,6 +83,8 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.5e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.91 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.005 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
# Parameters related to the initial conditions
InitialConditions:
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index 8ba338df28bffa5775d509fd5422913a8870e1d3..729be0c4ad628d727168cefc5a8e26e02844cac7 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -91,6 +91,8 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.5e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.91 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.005 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
# Parameters related to the initial conditions
InitialConditions:
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index e59d92389a6bfd526495ac70a008492b4d8ba8f2..e0e847d9c8e131a5211eb8fe7cb3941127ea3eb7 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -82,6 +82,8 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.5e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.91 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.005 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
# Parameters related to the initial conditions
InitialConditions:
diff --git a/examples/EAGLE_low_z/EAGLE_6/check_fof.py b/examples/EAGLE_low_z/EAGLE_6/check_fof.py
new file mode 100644
index 0000000000000000000000000000000000000000..4859af47497cebd50c46af3c8e441ff5f5af9829
--- /dev/null
+++ b/examples/EAGLE_low_z/EAGLE_6/check_fof.py
@@ -0,0 +1,472 @@
+import numpy as np
+import h5py as h5
+from tqdm import tqdm
+from numba import jit, prange
+
+snapname = "eagle_0000/eagle_0000.hdf5"
+fofname = "fof_output_0000/fof_output_0000.0.hdf5"
+# snapname = "eagle_0000.hdf5"
+# fofname = "fof_output_0000.hdf5"
+
+######################################################
+
+snap = h5.File(snapname, "r")
+
+nogrp_grp_id = int(snap["/Parameters"].attrs.get("FOF:group_id_default"))
+
+pos_gas = snap["/PartType0/Coordinates"][:, :]
+ids_gas = snap["/PartType0/ParticleIDs"][:]
+grp_gas = snap["/PartType0/FOFGroupIDs"][:]
+mass_gas = snap["/PartType0/Masses"][:]
+
+pos_DM = snap["/PartType1/Coordinates"][:, :]
+ids_DM = snap["/PartType1/ParticleIDs"][:]
+grp_DM = snap["/PartType1/FOFGroupIDs"][:]
+mass_DM = snap["/PartType1/Masses"][:]
+
+pos_star = snap["/PartType4/Coordinates"][:, :]
+ids_star = snap["/PartType4/ParticleIDs"][:]
+grp_star = snap["/PartType4/FOFGroupIDs"][:]
+mass_star = snap["/PartType4/Masses"][:]
+
+####################################################
+
+fof = h5.File(fofname, "r")
+num_files = fof["/Header/"].attrs["NumFilesPerSnapshot"][0]
+num_groups = fof["/Header/"].attrs["NumGroups_Total"][0]
+fof.close()
+
+fof_grp = np.zeros(num_groups, dtype=np.int32)
+fof_size = np.zeros(num_groups, dtype=np.int32)
+fof_mass = np.zeros(num_groups)
+
+# Read the distributed catalog
+offset = 0
+for i in range(num_files):
+
+ my_filename = fofname[:-6]
+ my_filename = my_filename + str(i) + ".hdf5"
+ fof = h5.File(my_filename, "r")
+
+ my_fof_grp = fof["/Groups/GroupIDs"][:]
+ my_fof_size = fof["/Groups/Sizes"][:]
+ my_fof_mass = fof["/Groups/Masses"][:]
+
+ num_this_file = fof["/Header"].attrs["NumGroups_ThisFile"][0]
+ fof.close()
+
+ fof_grp[offset : offset + num_this_file] = my_fof_grp
+ fof_size[offset : offset + num_this_file] = my_fof_size
+ fof_mass[offset : offset + num_this_file] = my_fof_mass
+
+ offset += num_this_file
+
+####################################################
+
+boxsize = snap["/Header"].attrs.get("BoxSize")[0]
+N_DM = snap["/Header"].attrs.get("NumPart_ThisFile")[1]
+
+l = 0.2 * boxsize / float(N_DM) ** (1.0 / 3.0)
+
+print("Checking snapshot :", snapname)
+print("Checking catalogue:", fofname)
+print("L:", boxsize)
+print("N_DM:", N_DM)
+print("Linking length:", l)
+print("")
+
+####################################################
+
+
+@jit(nopython=True, parallel=True, fastmath=True)
+def nearest(dx, L=boxsize):
+ mask1 = dx > 0.5 * L
+ mask2 = dx < -0.5 * L
+ if np.sum(mask1):
+ dx[mask1] = dx[mask1] - L
+ if np.sum(mask2):
+ dx[mask2] = dx[mask2] + L
+ return dx
+
+
+####################################################
+
+# Verify the content of the catalog
+num_groups = np.size(fof_grp)
+print("Catalog has", num_groups, "groups")
+
+
+def check_fof_size(i):
+ my_grp = fof_grp[i]
+ my_size = fof_size[i]
+
+ mask_gas = grp_gas == my_grp
+ mask_DM = grp_DM == my_grp
+ mask_star = grp_star == my_grp
+
+ total = np.sum(mask_gas) + np.sum(mask_DM) + np.sum(mask_star)
+
+ if total != my_size:
+ print(
+ "Grp",
+ my_grp,
+ "has size=",
+ my_size,
+ "but",
+ total,
+ "particles in the snapshot",
+ )
+ exit()
+
+
+for i in tqdm(range(num_groups)):
+ check_fof_size(i)
+
+print("All group sizes match the particles")
+####################################################
+
+# Verify group masses
+num_groups = np.size(fof_grp)
+print("Catalog has", num_groups, "groups")
+
+
+def check_fof_masses(i):
+ my_grp = fof_grp[i]
+ my_mass = fof_mass[i]
+
+ mask_gas = grp_gas == my_grp
+ mask_DM = grp_DM == my_grp
+ mask_star = grp_star == my_grp
+
+ total = (
+ np.sum(mass_gas[mask_gas])
+ + np.sum(mass_DM[mask_DM])
+ + np.sum(mass_star[mask_star])
+ )
+
+ ratio = total / my_mass
+
+ if ratio > 1.0001 or ratio < 0.9999:
+ print(
+ "Grp",
+ my_grp,
+ "has mass=",
+ my_mass,
+ "but particles in the snapshot have mass",
+ total,
+ )
+ exit()
+
+
+for i in tqdm(range(num_groups)):
+ check_fof_masses(i)
+
+print("All group masses match the particles")
+####################################################
+
+# Test the stand-alone stars
+mask = grp_star == nogrp_grp_id
+num_stars = np.sum(mask)
+print("Found %d stars not in groups" % num_stars)
+my_pos_star = pos_star[mask, :]
+my_ids_star = ids_star[mask]
+my_grp_star = grp_star[mask]
+my_pos_DM = pos_DM[:, :]
+my_ids_DM = ids_DM[:]
+my_grp_DM = grp_DM[:]
+
+# @jit(nopython=True, parallel=True, fastmath=True)
+def check_stand_alone_star(i):
+ pos = my_pos_star[i, :]
+ grp = my_grp_star[i]
+
+ dx = pos[0] - my_pos_DM[:, 0]
+ dy = pos[1] - my_pos_DM[:, 1]
+ dz = pos[2] - my_pos_DM[:, 2]
+
+ dx = nearest(dx)
+ dy = nearest(dy)
+ dz = nearest(dz)
+
+ # Identify the nearest DM particle
+ r2 = dx ** 2 + dy ** 2 + dz ** 2
+ select = np.argmin(r2)
+
+ # If the nearest DM particle is in a group --> mistake
+ target_grp = my_grp_DM[select]
+ if target_grp != nogrp_grp_id and r2[select] < l * l:
+ print("Found a star without group whose nearest DM particle is in a group!")
+ print("Star: id=", my_ids_star[i], "pos=", pos, "grp=", grp)
+ print(
+ "DM: id=",
+ my_ids_DM[select],
+ "pos=",
+ my_pos_DM[select],
+ "grp=",
+ my_grp_DM[select],
+ )
+ print("r=", np.sqrt(r2[select]))
+ # exit()
+
+
+for i in tqdm(range(num_stars)):
+ check_stand_alone_star(i)
+
+print("All stand-alone stars OK!")
+
+####################################################
+
+# Test the stars in groups
+mask = grp_star != nogrp_grp_id
+num_stars = np.sum(mask)
+print("Found %d stars in groups" % num_stars)
+my_pos_star = pos_star[mask, :]
+my_ids_star = ids_star[mask]
+my_grp_star = grp_star[mask]
+my_pos_DM = pos_DM[:, :]
+my_ids_DM = ids_DM[:]
+my_grp_DM = grp_DM[:]
+
+
+@jit(nopython=True, parallel=True, fastmath=True)
+def test_stars_in_group(i):
+ pos = my_pos_star[i, :]
+ grp = my_grp_star[i]
+
+ dx = pos[0] - my_pos_DM[:, 0]
+ dy = pos[1] - my_pos_DM[:, 1]
+ dz = pos[2] - my_pos_DM[:, 2]
+
+ dx = nearest(dx)
+ dy = nearest(dy)
+ dz = nearest(dz)
+
+ # Identify the nearest DM particle
+ r2 = dx ** 2 + dy ** 2 + dz ** 2
+ select = np.argmin(r2)
+
+ # If the nearest DM particle is not in the same group --> mistake
+ target_grp = my_grp_DM[select]
+ if target_grp != grp and r2[select] < l * l:
+ print(
+ "Found a star in a group whose nearest DM particle is in a different group!"
+ )
+ print("Star: id=", my_ids_star[i], "pos=", pos, "grp=", grp)
+ print(
+ "DM: id=", my_ids_DM[select], "pos=", my_pos_DM[select], "grp=", target_grp
+ )
+ print("r=", np.sqrt(r2[select]))
+ # exit()
+
+
+for i in tqdm(range(num_stars)):
+ test_stars_in_group(i)
+
+print("All stars in groups OK!")
+
+####################################################
+
+# Test the stand-alone gas
+mask = grp_gas == nogrp_grp_id
+num_gas = np.sum(mask)
+print("Found %d gas not in groups" % num_gas)
+my_pos_gas = pos_gas[mask, :]
+my_ids_gas = ids_gas[mask]
+my_grp_gas = grp_gas[mask]
+my_pos_DM = pos_DM[:, :]
+my_ids_DM = ids_DM[:]
+my_grp_DM = grp_DM[:]
+
+
+@jit(nopython=True, parallel=True, fastmath=True)
+def test_stand_alone_gas(i):
+ pos = my_pos_gas[i, :]
+ grp = my_grp_gas[i]
+
+ dx = pos[0] - my_pos_DM[:, 0]
+ dy = pos[1] - my_pos_DM[:, 1]
+ dz = pos[2] - my_pos_DM[:, 2]
+
+ dx = nearest(dx)
+ dy = nearest(dy)
+ dz = nearest(dz)
+
+ # Identify the nearest DM particle
+ r2 = dx ** 2 + dy ** 2 + dz ** 2
+ select = np.argmin(r2)
+
+ # If the nearest DM particle is in a group --> mistake
+ target_grp = my_grp_DM[select]
+ if target_grp != nogrp_grp_id and r2[select] < l * l:
+ print("Found a gas without group whose nearest DM particle is in a group!")
+ print("Gas: id=", my_ids_gas[i], "pos=", pos, "grp=", grp)
+ print(
+ "DM: id=",
+ my_ids_DM[select],
+ "pos=",
+ my_pos_DM[select],
+ "grp=",
+ my_grp_DM[select],
+ )
+ print("r=", np.sqrt(r2[select]))
+ # exit()
+
+
+for i in tqdm(range(num_gas)):
+ test_stand_alone_gas(i)
+
+print("All stand-alone gas OK!")
+
+####################################################
+
+# Test the gas in groups
+mask = grp_gas != nogrp_grp_id
+num_gas = np.sum(mask)
+print("Found %d gas in groups" % num_gas)
+my_pos_gas = pos_gas[mask, :]
+my_ids_gas = ids_gas[mask]
+my_grp_gas = grp_gas[mask]
+my_pos_DM = pos_DM[:, :]
+my_ids_DM = ids_DM[:]
+my_grp_DM = grp_DM[:]
+
+
+@jit(nopython=True, parallel=True, fastmath=True)
+def test_gas_in_groups(i):
+ pos = my_pos_gas[i, :]
+ grp = my_grp_gas[i]
+
+ dx = pos[0] - my_pos_DM[:, 0]
+ dy = pos[1] - my_pos_DM[:, 1]
+ dz = pos[2] - my_pos_DM[:, 2]
+
+ dx = nearest(dx)
+ dy = nearest(dy)
+ dz = nearest(dz)
+
+ # Identify the nearest DM particle
+ r2 = dx ** 2 + dy ** 2 + dz ** 2
+ select = np.argmin(r2)
+
+ # If the nearest DM particle is not in the same group --> mistake
+ target_grp = my_grp_DM[select]
+ if target_grp != grp and r2[select] < l * l:
+ print(
+ "Found a gas in a group whose nearest DM particle is in a different group!"
+ )
+ print("Gas: id=", my_ids_gas[i], "pos=", pos, "grp=", grp)
+ print(
+ "DM: id=", my_ids_DM[select], "pos=", my_pos_DM[select], "grp=", target_grp
+ )
+ print("r=", np.sqrt(r2[select]))
+ # exit()
+
+
+for i in tqdm(range(num_gas)):
+ test_gas_in_groups(i)
+
+print("All gas in groups OK!")
+
+####################################################
+
+# Test the stand-alone DM
+mask = grp_DM == nogrp_grp_id
+num_DM = np.sum(mask)
+print("Found %d DM not in groups" % num_DM)
+my_pos_DM = pos_DM[mask, :]
+my_ids_DM = ids_DM[mask]
+my_grp_DM = grp_DM[mask]
+
+
+@jit(nopython=True, parallel=True, fastmath=True)
+def test_stand_alone_DM(i):
+ pos = my_pos_DM[i, :]
+ grp = my_grp_DM[i]
+
+ dx = pos[0] - pos_DM[:, 0]
+ dy = pos[1] - pos_DM[:, 1]
+ dz = pos[2] - pos_DM[:, 2]
+
+ dx = nearest(dx)
+ dy = nearest(dy)
+ dz = nearest(dz)
+
+ # Identify the nearest DM particle
+ r2 = dx ** 2 + dy ** 2 + dz ** 2
+ mask = np.logical_and(r2 < l * l, r2 > 0.0)
+
+ # If the nearest DM particle is in a group --> mistake
+ if not np.all(grp_DM[mask] == nogrp_grp_id):
+ print("Found a DM without group with some DM particle within l in a group!")
+ print("DM: id=", my_ids_DM[i], "pos=", pos, "grp=", grp)
+ for j in range(np.sum(mask)):
+ if grp_DM[mask][j] != nogrp_grp_id:
+ print(
+ "Other: id=",
+ ids_DM[mask][j],
+ "pos=",
+ pos_DM[mask, :][j, :],
+ "grp=",
+ grp_DM[mask][j],
+ "r=",
+ np.sqrt(r2[mask][j]),
+ )
+
+
+for i in tqdm(range(num_DM)):
+ test_stand_alone_DM(i)
+
+print("All stand-alone DM OK!")
+
+####################################################
+
+# Test the DM in groups
+mask = grp_DM != nogrp_grp_id
+num_DM = np.sum(mask)
+print("Found %d DM in groups" % num_DM)
+my_pos_DM = pos_DM[mask, :]
+my_ids_DM = ids_DM[mask]
+my_grp_DM = grp_DM[mask]
+
+
+@jit(nopython=True, parallel=True, fastmath=True)
+def test_DM_in_groups(i):
+ pos = my_pos_DM[i, :]
+ grp = my_grp_DM[i]
+
+ dx = pos[0] - pos_DM[:, 0]
+ dy = pos[1] - pos_DM[:, 1]
+ dz = pos[2] - pos_DM[:, 2]
+
+ dx = nearest(dx)
+ dy = nearest(dy)
+ dz = nearest(dz)
+
+ # Identify the nearest DM particle
+ r2 = dx ** 2 + dy ** 2 + dz ** 2
+ mask = r2 < l * l
+
+ # If the nearest DM particle is not in the same group --> mistake
+ if not np.all(grp_DM[mask] == grp):
+ print(
+ "Found a DM in a group whose DM particles within l are in a different group!"
+ )
+ print("DM: id=", my_ids_DM[i], "pos=", pos, "grp=", grp)
+ for j in range(np.sum(mask)):
+ if grp_DM[mask][j] != grp:
+ print(
+ "Other: id=",
+ ids_DM[mask][j],
+ "pos=",
+ pos_DM[mask, :][j, :],
+ "grp=",
+ grp_DM[mask][j],
+ "r=",
+ np.sqrt(r2[mask][j]),
+ )
+
+
+for i in tqdm(range(num_DM)):
+ test_DM_in_groups(i)
+
+print("All DM in groups OK!")
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index a77e1bb827e8fbe9ad02e92e52bcfda1a25df28d..cf111232f1c20037dad7322e5ce64fad0069fb76 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -93,7 +93,9 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.5e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.91 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.005 # Scale-factor ratio between consecutive FoF black hole seeding calls.
-
+ linking_types: [0, 1, 0, 0, 0, 0, 0]
+ attaching_types: [1, 0, 0, 0, 1, 1, 0]
+
# Parameters related to the initial conditions
InitialConditions:
file_name: ./EAGLE_ICs_6.hdf5 # The file to read
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_lightcone/small_cosmo_volume.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_lightcone/small_cosmo_volume.yml
index 5bd006dc9739f2d21c580fa7182dc93929b2fb8d..31e662882e219c8b37bd2906995bde3eed91c3f2 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_lightcone/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_lightcone/small_cosmo_volume.yml
@@ -69,6 +69,8 @@ FOF:
black_hole_seed_halo_mass_Msun: 1.0e10 # Minimal halo mass in which to seed a black hole (in solar masses).
scale_factor_first: 0.05 # Scale-factor of first FoF black hole seeding calls.
delta_time: 1.00751 # Scale-factor ratio between consecutive FoF black hole seeding calls.
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
Scheduler:
max_top_level_cells: 8
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 0947c0da43928544de0f2a5f1f973748d8dec09e..2e16f8d7b938f79cca4fdca3a985f06f48a0b415 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -121,6 +121,8 @@ FOF:
group_id_offset: 1 # (Optional) Sets the offset of group ID labeling. Defaults to 1 if unspecified.
output_list_on: 0 # (Optional) Enable the output list
output_list: ./output_list_fof.txt # (Optional) File containing the output times (see documentation in "Parameter File" section)
+ linking_types: [0, 1, 0, 0, 0, 0, 0] # Use DM as the primary FOF linking type
+ attaching_types: [1, 0, 0, 0, 1, 1, 0] # Use gas, stars and black holes as FOF attachable types
# Parameters for the task scheduling
Scheduler:
diff --git a/src/common_io.c b/src/common_io.c
index 145eba9015802b8a5ecc871c16ddfd02bff37559..8d3941a246b35d4f490d12af1817bf6bb5ea41ab 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -507,10 +507,12 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str) {
* @param e The #engine containing the meta-data.
* @param internal_units The system of units used internally.
* @param snapshot_units The system of units used in snapshots.
+ * @param fof Is this a FOF output? If so don't write subgrid info.
*/
void io_write_meta_data(hid_t h_file, const struct engine* e,
const struct unit_system* internal_units,
- const struct unit_system* snapshot_units) {
+ const struct unit_system* snapshot_units,
+ const int fof) {
hid_t h_grp;
@@ -523,54 +525,57 @@ void io_write_meta_data(hid_t h_file, const struct engine* e,
/* Print the physical constants */
phys_const_print_snapshot(h_file, e->physical_constants);
- /* Print the SPH parameters */
- if (e->policy & engine_policy_hydro) {
- h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
+ if (!fof) {
+
+ /* Print the SPH parameters */
+ if (e->policy & engine_policy_hydro) {
+ h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
+ H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating SPH group");
+ hydro_props_print_snapshot(h_grp, e->hydro_properties);
+ hydro_write_flavour(h_grp);
+ mhd_write_flavour(h_grp);
+ H5Gclose(h_grp);
+ }
+
+ /* Print the subgrid parameters */
+ h_grp = H5Gcreate(h_file, "/SubgridScheme", H5P_DEFAULT, H5P_DEFAULT,
H5P_DEFAULT);
- if (h_grp < 0) error("Error while creating SPH group");
- hydro_props_print_snapshot(h_grp, e->hydro_properties);
- hydro_write_flavour(h_grp);
- mhd_write_flavour(h_grp);
+ if (h_grp < 0) error("Error while creating subgrid group");
+ hid_t h_grp_columns =
+ H5Gcreate(h_grp, "NamedColumns", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp_columns < 0) error("Error while creating named columns group");
+ entropy_floor_write_flavour(h_grp);
+ extra_io_write_flavour(h_grp, h_grp_columns);
+ cooling_write_flavour(h_grp, h_grp_columns, e->cooling_func);
+ chemistry_write_flavour(h_grp, h_grp_columns, e);
+ tracers_write_flavour(h_grp);
+ feedback_write_flavour(e->feedback_props, h_grp);
+ rt_write_flavour(h_grp, h_grp_columns, e, internal_units, snapshot_units,
+ e->rt_props);
H5Gclose(h_grp);
- }
- /* Print the subgrid parameters */
- h_grp = H5Gcreate(h_file, "/SubgridScheme", H5P_DEFAULT, H5P_DEFAULT,
- H5P_DEFAULT);
- if (h_grp < 0) error("Error while creating subgrid group");
- hid_t h_grp_columns =
- H5Gcreate(h_grp, "NamedColumns", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if (h_grp_columns < 0) error("Error while creating named columns group");
- entropy_floor_write_flavour(h_grp);
- extra_io_write_flavour(h_grp, h_grp_columns);
- cooling_write_flavour(h_grp, h_grp_columns, e->cooling_func);
- chemistry_write_flavour(h_grp, h_grp_columns, e);
- tracers_write_flavour(h_grp);
- feedback_write_flavour(e->feedback_props, h_grp);
- rt_write_flavour(h_grp, h_grp_columns, e, internal_units, snapshot_units,
- e->rt_props);
- H5Gclose(h_grp);
+ /* Print the gravity parameters */
+ if (e->policy & engine_policy_self_gravity) {
+ h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
+ H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating gravity group");
+ gravity_props_print_snapshot(h_grp, e->gravity_properties);
+ H5Gclose(h_grp);
+ }
- /* Print the gravity parameters */
- if (e->policy & engine_policy_self_gravity) {
- h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
- H5P_DEFAULT);
- if (h_grp < 0) error("Error while creating gravity group");
- gravity_props_print_snapshot(h_grp, e->gravity_properties);
- H5Gclose(h_grp);
- }
+ /* Print the stellar parameters */
+ if (e->policy & engine_policy_stars) {
+ h_grp = H5Gcreate(h_file, "/StarsScheme", H5P_DEFAULT, H5P_DEFAULT,
+ H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating stars group");
+ stars_props_print_snapshot(h_grp, h_grp_columns, e->stars_properties);
+ H5Gclose(h_grp);
+ }
- /* Print the stellar parameters */
- if (e->policy & engine_policy_stars) {
- h_grp = H5Gcreate(h_file, "/StarsScheme", H5P_DEFAULT, H5P_DEFAULT,
- H5P_DEFAULT);
- if (h_grp < 0) error("Error while creating stars group");
- stars_props_print_snapshot(h_grp, h_grp_columns, e->stars_properties);
- H5Gclose(h_grp);
+ H5Gclose(h_grp_columns);
}
- H5Gclose(h_grp_columns);
-
/* Print the cosmological model */
h_grp =
H5Gcreate(h_file, "/Cosmology", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
diff --git a/src/common_io.h b/src/common_io.h
index d7334a83fddf90aca12940cfe3290955fe65f10b..aced159e04e1bea23c5be84796efda9134896cc7 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -101,7 +101,8 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str);
void io_write_meta_data(hid_t h_file, const struct engine* e,
const struct unit_system* internal_units,
- const struct unit_system* snapshot_units);
+ const struct unit_system* snapshot_units,
+ const int fof);
void io_write_code_description(hid_t h_file);
void io_write_engine_policy(hid_t h_file, const struct engine* e);
diff --git a/src/distributed_io.c b/src/distributed_io.c
index b60958a3bea25239acf0c1d031b762efa66cef58..95cfa297f472ef4a450061c7bd24dd23d62374e9 100644
--- a/src/distributed_io.c
+++ b/src/distributed_io.c
@@ -452,6 +452,7 @@ void write_array_virtual(struct engine* e, hid_t grp, const char* fileName_base,
* @param numFields The number of fields to write for each particle type.
* @param internal_units The #unit_system used internally.
* @param snapshot_units The #unit_system used in the snapshots.
+ * @param fof Is this a snapshot related to a stand-alone FOF call?
* @param subsample_any Are any fields being subsampled?
* @param subsample_fraction The subsampling fraction of each particle type.
*/
@@ -463,7 +464,7 @@ void write_virtual_file(struct engine* e, const char* fileName_base,
const int numFields[swift_type_count],
char current_selection_name[FIELD_BUFFER_SIZE],
const struct unit_system* internal_units,
- const struct unit_system* snapshot_units,
+ const struct unit_system* snapshot_units, const int fof,
const int subsample_any,
const float subsample_fraction[swift_type_count]) {
@@ -604,7 +605,7 @@ void write_virtual_file(struct engine* e, const char* fileName_base,
ic_info_write_hdf5(e->ics_metadata, h_file);
/* Write all the meta-data */
- io_write_meta_data(h_file, e, internal_units, snapshot_units);
+ io_write_meta_data(h_file, e, internal_units, snapshot_units, fof);
/* Loop over all particle types */
for (int ptype = 0; ptype < swift_type_count; ptype++) {
@@ -738,6 +739,7 @@ void write_virtual_file(struct engine* e, const char* fileName_base,
* @param e The engine containing all the system.
* @param internal_units The #unit_system used internally
* @param snapshot_units The #unit_system used in the snapshots
+ * @param fof Is this a snapshot related to a stand-alone FOF call?
* @param mpi_rank The rank number of the calling MPI rank.
* @param mpi_size the number of MPI ranks.
* @param comm The communicator used by the MPI ranks.
@@ -751,8 +753,9 @@ void write_virtual_file(struct engine* e, const char* fileName_base,
void write_output_distributed(struct engine* e,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units,
- const int mpi_rank, const int mpi_size,
- MPI_Comm comm, MPI_Info info) {
+ const int fof, const int mpi_rank,
+ const int mpi_size, MPI_Comm comm,
+ MPI_Info info) {
hid_t h_file = 0, h_grp = 0;
int numFiles = mpi_size;
@@ -1097,7 +1100,7 @@ void write_output_distributed(struct engine* e,
ic_info_write_hdf5(e->ics_metadata, h_file);
/* Write all the meta-data */
- io_write_meta_data(h_file, e, internal_units, snapshot_units);
+ io_write_meta_data(h_file, e, internal_units, snapshot_units, fof);
/* Now write the top-level cell structure
* We use a global offset of 0 here. This means that the cells will write
@@ -1467,7 +1470,7 @@ void write_output_distributed(struct engine* e,
if (mpi_rank == 0)
write_virtual_file(e, fileName_base, xmfFileName, N_total, N_counts,
mpi_size, to_write, numFields, current_selection_name,
- internal_units, snapshot_units, subsample_any,
+ internal_units, snapshot_units, fof, subsample_any,
subsample_fraction);
/* Make sure nobody is allowed to progress until rank 0 is done. */
diff --git a/src/distributed_io.h b/src/distributed_io.h
index 0a73129ab7ab6c5a4502b4dcfeab9b72b47558f8..0fbb66cfe80e5a8fccf876a95983844334e6c9d0 100644
--- a/src/distributed_io.h
+++ b/src/distributed_io.h
@@ -35,8 +35,8 @@ struct unit_system;
void write_output_distributed(struct engine* e,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units,
- int mpi_rank, int mpi_size, MPI_Comm comm,
- MPI_Info info);
+ const int fof, int mpi_rank, int mpi_size,
+ MPI_Comm comm, MPI_Info info);
#endif /* HAVE_HDF5 && WITH_MPI */
diff --git a/src/engine.c b/src/engine.c
index b545fd66d1c2b1f92260c1010ad1a6119c5c6c48..bd71f4506f979dcd07f33ba4f45195f0e875bd81 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -782,12 +782,6 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) {
error(
"Not running with hydro but about to receive gas particles in "
"proxies!");
- if (!with_stars && count_sparts_in)
- error("Not running with stars but about to receive stars in proxies!");
- if (!with_black_holes && count_bparts_in)
- error(
- "Not running with black holes but about to receive black holes in "
- "proxies!");
if (e->verbose)
message("Counting number of foreign particles took %.3f %s.",
@@ -3628,7 +3622,7 @@ void engine_clean(struct engine *e, const int fof, const int restart) {
stats_free_mpi_type();
proxy_free_mpi_type();
task_free_mpi_comms();
- mpicollect_free_MPI_type();
+ if (!fof) mpicollect_free_MPI_type();
#endif
/* Close files */
diff --git a/src/engine.h b/src/engine.h
index 2760fc4c236c7e5ec58f2957e089dd3fde85acb9..c777508950068b1db069c2168f5a2d66134a75c8 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -696,7 +696,7 @@ void engine_io(struct engine *e);
void engine_io_check_snapshot_triggers(struct engine *e);
void engine_collect_end_of_step(struct engine *e, int apply);
void engine_collect_end_of_sub_cycle(struct engine *e);
-void engine_dump_snapshot(struct engine *e);
+void engine_dump_snapshot(struct engine *e, const int fof);
void engine_run_on_dump(struct engine *e);
void engine_init_output_lists(struct engine *e, struct swift_params *params,
const struct output_options *output_options);
@@ -760,6 +760,7 @@ void engine_fof(struct engine *e, const int dump_results,
const int dump_debug_results, const int seed_black_holes,
const int foreign_buffers_allocated);
void engine_activate_gpart_comms(struct engine *e);
+void engine_activate_fof_attach_tasks(struct engine *e);
/* Function prototypes, engine_maketasks.c. */
void engine_maketasks(struct engine *e);
diff --git a/src/engine_fof.c b/src/engine_fof.c
index 7cda72195d7cd10fee4b4102bb4c8569cbf9f87f..e98a50e263d08623d5a469ca4cdd03129aa058ef 100644
--- a/src/engine_fof.c
+++ b/src/engine_fof.c
@@ -68,7 +68,7 @@ void engine_activate_gpart_comms(struct engine *e) {
}
/**
- * @brief Activate all the FOF tasks.
+ * @brief Activate all the FOF linking tasks.
*
* Marks all the other task types to be skipped.
*
@@ -97,6 +97,37 @@ void engine_activate_fof_tasks(struct engine *e) {
clocks_getunit());
}
+/**
+ * @brief Activate all the FOF attaching tasks.
+ *
+ * Marks all the other task types to be skipped.
+ *
+ * @param e The #engine to act on.
+ */
+void engine_activate_fof_attach_tasks(struct engine *e) {
+
+ const ticks tic = getticks();
+
+ struct scheduler *s = &e->sched;
+ const int nr_tasks = s->nr_tasks;
+ struct task *tasks = s->tasks;
+
+ for (int k = 0; k < nr_tasks; k++) {
+
+ struct task *t = &tasks[k];
+
+ if (t->type == task_type_fof_attach_self ||
+ t->type == task_type_fof_attach_pair)
+ scheduler_activate(s, t);
+ else
+ t->skip = 1;
+ }
+
+ if (e->verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
+
/**
* @brief Run a FOF search.
*
@@ -123,15 +154,8 @@ void engine_fof(struct engine *e, const int dump_results,
#endif
}
- /* Compute number of DM particles */
- const long long total_nr_baryons =
- e->total_nr_parts + e->total_nr_sparts + e->total_nr_bparts;
- const long long total_nr_dmparts =
- e->total_nr_gparts - e->total_nr_DM_background_gparts -
- e->total_nr_neutrino_gparts - total_nr_baryons;
-
/* Initialise FOF parameters and allocate FOF arrays. */
- fof_allocate(e->s, total_nr_dmparts, e->fof_properties);
+ fof_allocate(e->s, e->fof_properties);
/* Make FOF tasks */
engine_make_fof_tasks(e);
@@ -142,14 +166,46 @@ void engine_fof(struct engine *e, const int dump_results,
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
- /* Perform local FOF tasks. */
+ /* Perform local FOF tasks for linkable particles. */
engine_launch(e, "fof");
- /* Perform FOF search over foreign particles and
- * find groups which require black hole seeding. */
- fof_search_tree(e->fof_properties, e->black_holes_properties,
- e->physical_constants, e->cosmology, e->s, dump_results,
- dump_debug_results, seed_black_holes);
+ /* Compute group sizes (only of local fragments with MPI) */
+ fof_compute_local_sizes(e->fof_properties, e->s);
+
+#ifdef WITH_MPI
+
+ /* Allocate buffers to receive the gpart fof information */
+ engine_allocate_foreign_particles(e, /*fof=*/1);
+
+ /* Compute the local<->foreign group links (nothing to do without MPI)*/
+ fof_search_foreign_cells(e->fof_properties, e->s);
+#endif
+
+ /* Compute the attachable->linkable links */
+ fof_link_attachable_particles(e->fof_properties, e->s);
+
+#ifdef WITH_MPI
+
+ /* Free the foreign particles */
+ space_free_foreign_parts(e->s, /*clear pointers=*/1);
+
+#endif
+
+ /* Finish the operations attaching the attachables to their groups */
+ fof_finalise_attachables(e->fof_properties, e->s);
+
+#ifdef WITH_MPI
+
+ /* Link the foreign fragments and finalise global group list (nothing to do
+ * without MPI) */
+ fof_link_foreign_fragments(e->fof_properties, e->s);
+#endif
+
+ /* Compute group properties and act on the results
+ * (seed BHs, dump catalogues..) */
+ fof_compute_group_props(e->fof_properties, e->black_holes_properties,
+ e->physical_constants, e->cosmology, e->s,
+ dump_results, dump_debug_results, seed_black_holes);
/* Reset flag. */
e->run_fof = 0;
diff --git a/src/engine_io.c b/src/engine_io.c
index 279e341153d85c0aa856b9ba6a59e719990a8044..a9291bbd618c15b1fcd334f1d2e75cbf8e870182 100644
--- a/src/engine_io.c
+++ b/src/engine_io.c
@@ -268,8 +268,9 @@ int engine_dump_restarts(struct engine *e, const int drifted_all,
* @brief Writes a snapshot with the current state of the engine
*
* @param e The #engine.
+ * @param fof Is this a stand-alone FOF call?
*/
-void engine_dump_snapshot(struct engine *e) {
+void engine_dump_snapshot(struct engine *e, const int fof) {
struct clocks_time time1, time2;
clocks_gettime(&time1);
@@ -326,20 +327,22 @@ void engine_dump_snapshot(struct engine *e) {
if (e->snapshot_distributed) {
- write_output_distributed(e, e->internal_units, e->snapshot_units, e->nodeID,
- e->nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
+ write_output_distributed(e, e->internal_units, e->snapshot_units, fof,
+ e->nodeID, e->nr_nodes, MPI_COMM_WORLD,
+ MPI_INFO_NULL);
} else {
#if defined(HAVE_PARALLEL_HDF5)
- write_output_parallel(e, e->internal_units, e->snapshot_units, e->nodeID,
- e->nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
+ write_output_parallel(e, e->internal_units, e->snapshot_units, fof,
+ e->nodeID, e->nr_nodes, MPI_COMM_WORLD,
+ MPI_INFO_NULL);
#else
- write_output_serial(e, e->internal_units, e->snapshot_units, e->nodeID,
+ write_output_serial(e, e->internal_units, e->snapshot_units, fof, e->nodeID,
e->nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
#endif
}
#else
- write_output_single(e, e->internal_units, e->snapshot_units);
+ write_output_single(e, e->internal_units, e->snapshot_units, fof);
#endif
#endif
@@ -541,7 +544,7 @@ void engine_io(struct engine *e) {
}
/* Dump... */
- engine_dump_snapshot(e);
+ engine_dump_snapshot(e, /*fof=*/0);
/* Free the memory allocated for VELOCIraptor i/o. */
if (with_stf && e->snapshot_invoke_stf && e->s->gpart_group_data) {
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index e205cf1f1a74371d7b117868095328b071bbb51f..376c6966a6aad7fe4fa89ba78157619f1a488c5d 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -4490,11 +4490,12 @@ void engine_make_fofloop_tasks_mapper(void *map_data, int num_elements,
if (ci->grav.count == 0) continue;
/* If the cells is local build a self-interaction */
- if (ci->nodeID == nodeID)
+ if (ci->nodeID == nodeID) {
scheduler_addtask(sched, task_type_fof_self, task_subtype_none, 0, 0, ci,
NULL);
- else
- continue;
+ scheduler_addtask(sched, task_type_fof_attach_self, task_subtype_none, 0,
+ 0, ci, NULL);
+ }
/* Now loop over all the neighbours of this cell */
for (int ii = -1; ii < 2; ii++) {
@@ -4514,13 +4515,19 @@ void engine_make_fofloop_tasks_mapper(void *map_data, int num_elements,
const int cjd = cell_getid(cdim, iii, jjj, kkk);
struct cell *cj = &cells[cjd];
- /* Is that neighbour local and does it have particles ? */
- if (cid >= cjd || cj->grav.count == 0 || (ci->nodeID != cj->nodeID))
- continue;
+ /* Does that neighbour have particles ? */
+ if (cid >= cjd || cj->grav.count == 0) continue;
- /* Construct the pair task */
- scheduler_addtask(sched, task_type_fof_pair, task_subtype_none, 0, 0,
- ci, cj);
+ /* Construct the pair search task only for fully local pairs */
+ if (ci->nodeID == nodeID && cj->nodeID == nodeID)
+ scheduler_addtask(sched, task_type_fof_pair, task_subtype_none, 0,
+ 0, ci, cj);
+
+ /* Construct the pair search task for pairs overlapping with the node
+ */
+ if (ci->nodeID == nodeID || cj->nodeID == nodeID)
+ scheduler_addtask(sched, task_type_fof_attach_pair,
+ task_subtype_none, 0, 0, ci, cj);
}
}
}
diff --git a/src/fof.c b/src/fof.c
index 50c2f726cb944a6677097e96c72c9b4b525592c2..4c9a37c514ac84a8053786fbd53fd6e6a15d4e13 100644
--- a/src/fof.c
+++ b/src/fof.c
@@ -55,6 +55,15 @@
#define UNION_BY_SIZE_OVER_MPI (1)
#define FOF_COMPRESS_PATHS_MIN_LENGTH (2)
+/* The FoF policy we are running */
+int current_fof_linking_type;
+
+/* The FoF policy for particles attached to the main type */
+int current_fof_attach_type;
+
+/* The FoF policy for particles ignored altogether */
+int current_fof_ignore_type;
+
/* Are we timing calculating group properties in the FOF? */
//#define WITHOUT_GROUP_PROPS
@@ -148,6 +157,62 @@ void fof_init(struct fof_props *props, struct swift_params *params,
props->seed_halo_mass *= phys_const->const_solar_mass;
}
+ /* Read what particle types we want to run FOF on */
+ parser_get_param_int_array(params, "FOF:linking_types", swift_type_count,
+ props->fof_linking_types);
+
+ /* Read what particle types we want to attach to FOF groups */
+ parser_get_param_int_array(params, "FOF:attaching_types", swift_type_count,
+ props->fof_attach_types);
+
+ /* Check that there is something to do */
+ int sum = 0;
+ for (int i = 0; i < swift_type_count; ++i)
+ if (props->fof_linking_types[i]) sum++;
+ if (sum == 0) error("FOF must run on at least one type of particles!");
+
+ for (int i = 0; i < swift_type_count; ++i)
+ if (props->fof_linking_types[i] && props->fof_attach_types[i])
+ error("FOF can't use a type (%s) as both linking and attaching type!",
+ part_type_names[i]);
+
+ /* Initialize the FoF linking mode */
+ current_fof_linking_type = 0;
+ for (int i = 0; i < swift_type_count; ++i)
+ if (props->fof_linking_types[i]) {
+ current_fof_linking_type |= (1 << (i + 1));
+ }
+
+ /* Initialize the FoF attaching mode */
+ current_fof_attach_type = 0;
+ for (int i = 0; i < swift_type_count; ++i)
+ if (props->fof_attach_types[i]) {
+ current_fof_attach_type |= (1 << (i + 1));
+ }
+
+ /* Construct the combined mask of ignored particles */
+ current_fof_ignore_type =
+ ~(current_fof_linking_type | current_fof_attach_type);
+
+ /* Report what we do */
+ if (engine_rank == 0) {
+ printf("FOF using the following types for linking:");
+ for (int i = 0; i < swift_type_count; ++i)
+ if (props->fof_linking_types[i]) printf("'%s' ", part_type_names[i]);
+ printf("\n");
+
+ printf("FOF using the following types for attaching:");
+ for (int i = 0; i < swift_type_count; ++i)
+ if (props->fof_attach_types[i]) printf("'%s' ", part_type_names[i]);
+ printf("\n");
+
+ printf("FOF ignoring the following types:");
+ for (int i = 0; i < swift_type_count; ++i)
+ if (current_fof_ignore_type & (1 << (i + 1)))
+ printf("'%s' ", part_type_names[i]);
+ printf("\n");
+ }
+
#if defined(WITH_MPI) && defined(UNION_BY_SIZE_OVER_MPI)
if (engine_rank == 0)
message(
@@ -210,6 +275,25 @@ void fof_set_initial_group_index_mapper(void *map_data, int num_elements,
}
}
+/**
+ * @brief Mapper function to set the initial attach group indices.
+ *
+ * @param map_data The array of attach group indices.
+ * @param num_elements Chunk size.
+ * @param extra_data Pointer to first group index.
+ */
+void fof_set_initial_attach_index_mapper(void *map_data, int num_elements,
+ void *extra_data) {
+ size_t *attach_index = (size_t *)map_data;
+ size_t *attach_index_start = (size_t *)extra_data;
+
+ const ptrdiff_t offset = attach_index - attach_index_start;
+
+ for (int i = 0; i < num_elements; ++i) {
+ attach_index[i] = i + offset;
+ }
+}
+
/**
* @brief Mapper function to set the initial group sizes.
*
@@ -226,6 +310,22 @@ void fof_set_initial_group_size_mapper(void *map_data, int num_elements,
}
}
+/**
+ * @brief Mapper function to set the initial distances.
+ *
+ * @param map_data The array of distance.
+ * @param num_elements Chunk size.
+ * @param extra_data N/A.
+ */
+void fof_set_initial_part_distances_mapper(void *map_data, int num_elements,
+ void *extra_data) {
+
+ float *distance = (float *)map_data;
+ for (int i = 0; i < num_elements; ++i) {
+ distance[i] = FLT_MAX;
+ }
+}
+
/**
* @brief Mapper function to set the initial group IDs.
*
@@ -249,30 +349,38 @@ void fof_set_initial_group_id_mapper(void *map_data, int num_elements,
* @brief Allocate the memory and initialise the arrays for a FOF calculation.
*
* @param s The #space to act on.
- * @param total_nr_DM_particles The total number of DM particles in the
- * simulation.
* @param props The properties of the FOF structure.
*/
-void fof_allocate(const struct space *s, const long long total_nr_DM_particles,
- struct fof_props *props) {
+void fof_allocate(const struct space *s, struct fof_props *props) {
const int verbose = s->e->verbose;
const ticks total_tic = getticks();
/* Start by computing the mean inter DM particle separation */
- /* Collect the mass of the first non-background gpart */
+ /* Collect the mean mass of the non-background gpart */
double high_res_DM_mass = 0.;
+ long long num_high_res_DM = 0;
for (size_t i = 0; i < s->nr_gparts; ++i) {
const struct gpart *gp = &s->gparts[i];
if (gp->type == swift_type_dark_matter &&
gp->time_bin != time_bin_inhibited &&
gp->time_bin != time_bin_not_created) {
- high_res_DM_mass = gp->mass;
- break;
+ high_res_DM_mass += gp->mass;
+ num_high_res_DM++;
}
}
+#ifdef WITH_MPI
+ /* Gather the information from all ranks */
+ MPI_Allreduce(MPI_IN_PLACE, &high_res_DM_mass, 1, MPI_DOUBLE, MPI_SUM,
+ MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, &num_high_res_DM, 1, MPI_LONG_LONG, MPI_SUM,
+ MPI_COMM_WORLD);
+#endif
+
+ high_res_DM_mass /= (double)num_high_res_DM;
+
/* Are we using the aboslute value or the one derived from the mean
inter-particle sepration? */
if (props->l_x_absolute != -1.) {
@@ -334,6 +442,24 @@ void fof_allocate(const struct space *s, const long long total_nr_DM_particles,
s->nr_gparts * sizeof(size_t)) != 0)
error("Failed to allocate list of particle group indices for FOF search.");
+ /* Allocate and initialise a group index array for attachables. */
+ if (swift_memalign("fof_attach_index", (void **)&props->attach_index, 64,
+ s->nr_gparts * sizeof(size_t)) != 0)
+ error(
+ "Failed to allocate list of particle distances array for FOF search.");
+
+ /* Allocate and initialise a group index array for attachables. */
+ if (swift_memalign("fof_found_attach", (void **)&props->found_attachable_link,
+ 64, s->nr_gparts * sizeof(char)) != 0)
+ error(
+ "Failed to allocate list of particle distances array for FOF search.");
+
+ /* Allocate and initialise the closest distance array. */
+ if (swift_memalign("fof_distance", (void **)&props->distance_to_link, 64,
+ s->nr_gparts * sizeof(float)) != 0)
+ error(
+ "Failed to allocate list of particle distances array for FOF search.");
+
/* Allocate and initialise a group size array. */
if (swift_memalign("fof_group_size", (void **)&props->group_size, 64,
s->nr_gparts * sizeof(size_t)) != 0)
@@ -352,6 +478,30 @@ void fof_allocate(const struct space *s, const long long total_nr_DM_particles,
tic = getticks();
+ /* Set initial attach index */
+ threadpool_map(&s->e->threadpool, fof_set_initial_attach_index_mapper,
+ props->attach_index, s->nr_gparts, sizeof(size_t),
+ threadpool_auto_chunk_size, props->attach_index);
+
+ bzero(props->found_attachable_link, s->nr_gparts * sizeof(char));
+
+ if (verbose)
+ message("Setting initial attach index took: %.3f %s.",
+ clocks_from_ticks(getticks() - tic), clocks_getunit());
+
+ tic = getticks();
+
+ /* Set initial distances */
+ threadpool_map(&s->e->threadpool, fof_set_initial_part_distances_mapper,
+ props->distance_to_link, s->nr_gparts, sizeof(float),
+ threadpool_auto_chunk_size, NULL);
+
+ if (verbose)
+ message("Setting initial distances took: %.3f %s.",
+ clocks_from_ticks(getticks() - tic), clocks_getunit());
+
+ tic = getticks();
+
/* Set initial group sizes */
threadpool_map(&s->e->threadpool, fof_set_initial_group_size_mapper,
props->group_size, s->nr_gparts, sizeof(size_t),
@@ -491,7 +641,6 @@ __attribute__((always_inline)) INLINE static size_t fof_find_global(
#endif /* WITH_MPI */
-#ifndef WITHOUT_GROUP_PROPS
/**
* @brief Finds the local root ID of the group a particle exists in
* when group_index contains globally unique identifiers -
@@ -526,7 +675,33 @@ __attribute__((always_inline)) INLINE static size_t fof_find_local(
return root;
#endif
}
-#endif /* #ifndef WITHOUT_GROUP_PROPS */
+
+/**
+ * @brief Returns whether a #gpart is of the 'attachable' kind.
+ */
+__attribute__((always_inline)) INLINE static int gpart_is_attachable(
+ const struct gpart *gp) {
+
+ return current_fof_attach_type & (1 << (gp->type + 1));
+}
+
+/**
+ * @brief Returns whether a #gpart is of the 'linkable' kind.
+ */
+__attribute__((always_inline)) INLINE static int gpart_is_linkable(
+ const struct gpart *gp) {
+
+ return current_fof_linking_type & (1 << (gp->type + 1));
+}
+
+/**
+ * @brief Returns whether a #gpart is to be ignored by FOF.
+ */
+__attribute__((always_inline)) INLINE static int gpart_is_ignorable(
+ const struct gpart *gp) {
+
+ return current_fof_ignore_type & (1 << (gp->type + 1));
+}
/**
* @brief Finds the local root ID of the group a particle exists in.
@@ -596,7 +771,8 @@ __attribute__((always_inline)) INLINE static int atomic_update_root(
* @param group_index The list of group roots.
*/
__attribute__((always_inline)) INLINE static void fof_union(
- size_t *root_i, const size_t root_j, size_t *group_index) {
+ size_t *restrict root_i, const size_t root_j,
+ size_t *restrict group_index) {
int result = 0;
@@ -789,34 +965,36 @@ __attribute__((always_inline)) INLINE static void fof_compute_send_recv_offsets(
*/
void fof_search_self_cell(const struct fof_props *props, const double l_x2,
const struct gpart *const space_gparts,
- struct cell *c) {
+ const struct cell *c) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->split) error("Performing the FOF search at a non-leaf level!");
#endif
const size_t count = c->grav.count;
- struct gpart *gparts = c->grav.parts;
+ const struct gpart *gparts = c->grav.parts;
/* Index of particles in the global group list */
- size_t *group_index = props->group_index;
+ size_t *const group_index = props->group_index;
/* Make a list of particle offsets into the global gparts array. */
size_t *const offset = group_index + (ptrdiff_t)(gparts - space_gparts);
+#ifdef SWIFT_DEBUG_CHECKS
if (c->nodeID != engine_rank)
error("Performing self FOF search on foreign cell.");
+#endif
/* Loop over particles and find which particles belong in the same group. */
for (size_t i = 0; i < count; i++) {
- struct gpart *pi = &gparts[i];
+ const struct gpart *pi = &gparts[i];
/* Ignore inhibited particles */
if (pi->time_bin >= time_bin_inhibited) continue;
- /* Ignore neutrinos */
- if (pi->type == swift_type_neutrino) continue;
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(pi)) continue;
#ifdef SWIFT_DEBUG_CHECKS
if (pi->ti_drift != ti_current)
@@ -830,31 +1008,40 @@ void fof_search_self_cell(const struct fof_props *props, const double l_x2,
/* Find the root of pi. */
size_t root_i = fof_find(offset[i], group_index);
+ /* Get the nature of the linking */
+ const int is_link_i = gpart_is_linkable(pi);
+
for (size_t j = i + 1; j < count; j++) {
- struct gpart *pj = &gparts[j];
+ const struct gpart *pj = &gparts[j];
/* Ignore inhibited particles */
if (pj->time_bin >= time_bin_inhibited) continue;
- /* Ignore neutrinos */
- if (pj->type == swift_type_neutrino) continue;
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(pj)) continue;
+
+ /* Get the nature of the linking */
+ const int is_link_j = gpart_is_linkable(pj);
+
+ /* Both particles must be of the linking kind */
+ if (!(is_link_i && is_link_j)) continue;
#ifdef SWIFT_DEBUG_CHECKS
if (pj->ti_drift != ti_current)
error("Running FOF on an un-drifted particle!");
#endif
- const double pjx = pj->x[0];
- const double pjy = pj->x[1];
- const double pjz = pj->x[2];
-
/* Find the root of pj. */
const size_t root_j = fof_find(offset[j], group_index);
/* Skip particles in the same group. */
if (root_i == root_j) continue;
+ const double pjx = pj->x[0];
+ const double pjy = pj->x[1];
+ const double pjz = pj->x[2];
+
/* Compute the pairwise distance */
float dx[3], r2 = 0.0f;
dx[0] = pix - pjx;
@@ -866,7 +1053,7 @@ void fof_search_self_cell(const struct fof_props *props, const double l_x2,
/* Hit or miss? */
if (r2 < l_x2) {
- /* Merge the groups */
+ /* Merge the groups` */
fof_union(&root_i, root_j, group_index);
}
}
@@ -887,15 +1074,16 @@ void fof_search_self_cell(const struct fof_props *props, const double l_x2,
void fof_search_pair_cells(const struct fof_props *props, const double dim[3],
const double l_x2, const int periodic,
const struct gpart *const space_gparts,
- struct cell *restrict ci, struct cell *restrict cj) {
+ const struct cell *restrict ci,
+ const struct cell *restrict cj) {
const size_t count_i = ci->grav.count;
const size_t count_j = cj->grav.count;
- struct gpart *gparts_i = ci->grav.parts;
- struct gpart *gparts_j = cj->grav.parts;
+ const struct gpart *gparts_i = ci->grav.parts;
+ const struct gpart *gparts_j = cj->grav.parts;
/* Index of particles in the global group list */
- size_t *group_index = props->group_index;
+ size_t *const group_index = props->group_index;
/* Make a list of particle offsets into the global gparts array. */
size_t *const offset_i = group_index + (ptrdiff_t)(gparts_i - space_gparts);
@@ -928,13 +1116,13 @@ void fof_search_pair_cells(const struct fof_props *props, const double dim[3],
/* Loop over particles and find which particles belong in the same group. */
for (size_t i = 0; i < count_i; i++) {
- struct gpart *restrict pi = &gparts_i[i];
+ const struct gpart *restrict pi = &gparts_i[i];
/* Ignore inhibited particles */
if (pi->time_bin >= time_bin_inhibited) continue;
- /* Ignore neutrinos */
- if (pi->type == swift_type_neutrino) continue;
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(pi)) continue;
#ifdef SWIFT_DEBUG_CHECKS
if (pi->ti_drift != ti_current)
@@ -948,15 +1136,24 @@ void fof_search_pair_cells(const struct fof_props *props, const double dim[3],
/* Find the root of pi. */
size_t root_i = fof_find(offset_i[i], group_index);
+ /* Get the nature of the linking */
+ const int is_link_i = gpart_is_linkable(pi);
+
for (size_t j = 0; j < count_j; j++) {
- struct gpart *restrict pj = &gparts_j[j];
+ const struct gpart *restrict pj = &gparts_j[j];
/* Ignore inhibited particles */
if (pj->time_bin >= time_bin_inhibited) continue;
- /* Ignore neutrinos */
- if (pj->type == swift_type_neutrino) continue;
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(pj)) continue;
+
+ /* Get the nature of the linking */
+ const int is_link_j = gpart_is_linkable(pj);
+
+ /* At least one of the particles has to be of linking type */
+ if (!(is_link_i && is_link_j)) continue;
#ifdef SWIFT_DEBUG_CHECKS
if (pj->ti_drift != ti_current)
@@ -973,8 +1170,8 @@ void fof_search_pair_cells(const struct fof_props *props, const double dim[3],
const double pjy = pj->x[1];
const double pjz = pj->x[2];
- /* Compute pairwise distance, remembering to account for boundary
- * conditions. */
+ /* Compute pairwise distance (periodic BCs were accounted
+ for by the shift vector) */
float dx[3], r2 = 0.0f;
dx[0] = pix - pjx;
dx[1] = piy - pjy;
@@ -992,6 +1189,48 @@ void fof_search_pair_cells(const struct fof_props *props, const double dim[3],
}
}
+#ifdef WITH_MPI
+
+/**
+ * @brief Add a local<->foreign pair in range to the list of links
+ *
+ * Possibly reallocates the local_group_links if we run out of space.
+ */
+static INLINE void add_foreign_link_to_list(
+ int *local_link_count, int *group_links_size, struct fof_mpi **group_links,
+ struct fof_mpi **local_group_links, const size_t root_i,
+ const size_t root_j, const size_t size_i, const size_t size_j) {
+
+ /* If the group_links array is not big enough re-allocate it. */
+ if (*local_link_count + 1 > *group_links_size) {
+
+ const int new_size = 2 * (*group_links_size);
+
+ *group_links_size = new_size;
+
+ (*group_links) = (struct fof_mpi *)realloc(
+ *group_links, new_size * sizeof(struct fof_mpi));
+
+ /* Reset the local pointer */
+ (*local_group_links) = *group_links;
+
+ message("Re-allocating local group links from %d to %d elements.",
+ *local_link_count, new_size);
+ }
+
+ /* Store the particle group properties for communication. */
+
+ (*local_group_links)[*local_link_count].group_i = root_i;
+ (*local_group_links)[*local_link_count].group_i_size = size_i;
+
+ (*local_group_links)[*local_link_count].group_j = root_j;
+ (*local_group_links)[*local_link_count].group_j_size = size_j;
+
+ (*local_link_count)++;
+}
+
+#endif
+
/* Perform a FOF search between a local and foreign cell using the Union-Find
* algorithm. Store any links found between particles.*/
void fof_search_pair_cells_foreign(
@@ -1008,15 +1247,16 @@ void fof_search_pair_cells_foreign(
const struct gpart *gparts_j = cj->grav.parts;
/* Get local pointers */
- size_t *group_index = props->group_index;
- size_t *group_size = props->group_size;
+ const size_t *restrict group_index = props->group_index;
+ const size_t *restrict group_size = props->group_size;
/* Values local to this function to avoid dereferencing */
struct fof_mpi *local_group_links = *group_links;
int local_link_count = *link_count;
/* Make a list of particle offsets into the global gparts array. */
- size_t *const offset_i = group_index + (ptrdiff_t)(gparts_i - space_gparts);
+ const size_t *const offset_i =
+ group_index + (ptrdiff_t)(gparts_i - space_gparts);
#ifdef SWIFT_DEBUG_CHECKS
@@ -1055,8 +1295,8 @@ void fof_search_pair_cells_foreign(
/* Ignore inhibited particles */
if (pi->time_bin >= time_bin_inhibited) continue;
- /* Ignore neutrinos */
- if (pi->type == swift_type_neutrino) continue;
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(pi)) continue;
#ifdef SWIFT_DEBUG_CHECKS
if (pi->ti_drift != ti_current)
@@ -1071,6 +1311,9 @@ void fof_search_pair_cells_foreign(
const size_t root_i =
fof_find_global(offset_i[i] - node_offset, group_index, nr_gparts);
+ /* Get the nature of the linking */
+ const int is_link_i = gpart_is_linkable(pi);
+
for (size_t j = 0; j < count_j; j++) {
const struct gpart *pj = &gparts_j[j];
@@ -1078,8 +1321,14 @@ void fof_search_pair_cells_foreign(
/* Ignore inhibited particles */
if (pj->time_bin >= time_bin_inhibited) continue;
- /* Ignore neutrinos */
- if (pj->type == swift_type_neutrino) continue;
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(pj)) continue;
+
+ /* Get the nature of the linking */
+ const int is_link_j = gpart_is_linkable(pj);
+
+ /* Only consider linkable<->linkable pairs */
+ if (!(is_link_i && is_link_j)) continue;
#ifdef SWIFT_DEBUG_CHECKS
if (pj->ti_drift != ti_current)
@@ -1090,8 +1339,8 @@ void fof_search_pair_cells_foreign(
const double pjy = pj->x[1];
const double pjz = pj->x[2];
- /* Compute pairwise distance, remembering to account for boundary
- * conditions. */
+ /* Compute pairwise distance (periodic BCs were accounted
+ for by the shift vector) */
float dx[3], r2 = 0.0f;
dx[0] = pix - pjx;
dx[1] = piy - pjy;
@@ -1102,47 +1351,19 @@ void fof_search_pair_cells_foreign(
/* Hit or miss? */
if (r2 < l_x2) {
- int found = 0;
-
/* Check that the links have not already been added to the list. */
for (int l = 0; l < local_link_count; l++) {
- if ((local_group_links)[l].group_i == root_i &&
- (local_group_links)[l].group_j == pj->fof_data.group_id) {
- found = 1;
- break;
+ if (local_group_links[l].group_i == root_i &&
+ local_group_links[l].group_j == pj->fof_data.group_id) {
+ continue;
}
}
- if (!found) {
-
- /* If the group_links array is not big enough re-allocate it. */
- if (local_link_count + 1 > *group_links_size) {
-
- const int new_size = 2 * (*group_links_size);
-
- *group_links_size = new_size;
-
- (*group_links) = (struct fof_mpi *)realloc(
- *group_links, new_size * sizeof(struct fof_mpi));
-
- /* Reset the local pointer */
- local_group_links = *group_links;
-
- message("Re-allocating local group links from %d to %d elements.",
- local_link_count, new_size);
- }
-
- /* Store the particle group properties for communication. */
- local_group_links[local_link_count].group_i = root_i;
- local_group_links[local_link_count].group_i_size =
- group_size[root_i - node_offset];
-
- local_group_links[local_link_count].group_j = pj->fof_data.group_id;
- local_group_links[local_link_count].group_j_size =
- pj->fof_data.group_size;
-
- local_link_count++;
- }
+ /* Add a possible link to the list */
+ add_foreign_link_to_list(
+ &local_link_count, group_links_size, group_links,
+ &local_group_links, root_i, pj->fof_data.group_id,
+ group_size[root_i - node_offset], pj->fof_data.group_size);
}
}
}
@@ -1225,17 +1446,546 @@ void rec_fof_search_pair_foreign(
int *restrict link_count, struct fof_mpi **group_links,
int *restrict group_links_size) {
-#ifdef SWIFT_DEBUG_CHECKS
- if (ci == cj) error("Pair FOF called on same cell!!!");
- if (ci->nodeID == cj->nodeID) error("Fully local pair!");
-#endif
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci == cj) error("Pair FOF called on same cell!!!");
+ if (ci->nodeID == cj->nodeID) error("Fully local pair!");
+#endif
+
+ /* Find the shortest distance between cells, remembering to account for
+ * boundary conditions. */
+ const double r2 = cell_min_dist(ci, cj, dim);
+
+ /* Return if cells are out of range of each other. */
+ if (r2 > search_r2) return;
+
+ /* Recurse on both cells if they are both split. */
+ if (ci->split && cj->split) {
+ for (int k = 0; k < 8; k++) {
+ if (ci->progeny[k] != NULL) {
+
+ for (int l = 0; l < 8; l++)
+ if (cj->progeny[l] != NULL)
+ rec_fof_search_pair_foreign(props, dim, search_r2, periodic,
+ space_gparts, nr_gparts, ci->progeny[k],
+ cj->progeny[l], link_count, group_links,
+ group_links_size);
+ }
+ }
+ } else if (ci->split) {
+
+ for (int k = 0; k < 8; k++) {
+ if (ci->progeny[k] != NULL)
+ rec_fof_search_pair_foreign(props, dim, search_r2, periodic,
+ space_gparts, nr_gparts, ci->progeny[k], cj,
+ link_count, group_links, group_links_size);
+ }
+ } else if (cj->split) {
+ for (int k = 0; k < 8; k++) {
+ if (cj->progeny[k] != NULL)
+ rec_fof_search_pair_foreign(props, dim, search_r2, periodic,
+ space_gparts, nr_gparts, ci, cj->progeny[k],
+ link_count, group_links, group_links_size);
+ }
+ } else {
+ /* Perform FOF search between pairs of cells that are within the linking
+ * length and not the same cell. */
+ fof_search_pair_cells_foreign(props, dim, search_r2, periodic, space_gparts,
+ nr_gparts, ci, cj, link_count, group_links,
+ group_links_size);
+ }
+}
+
+#endif /* WITH_MPI */
+
+/**
+ * @brief Recursively perform a union-find FOF on a cell.
+ *
+ * @param props The properties fof the FOF scheme.
+ * @param dim The dimension of the space.
+ * @param space_gparts The start of the #gpart array in the #space structure.
+ * @param search_r2 the square of the FOF linking length.
+ * @param periodic Are we using periodic BCs?
+ * @param c The #cell in which to perform FOF.
+ */
+void rec_fof_search_self(const struct fof_props *props, const double dim[3],
+ const double search_r2, const int periodic,
+ const struct gpart *const space_gparts,
+ struct cell *c) {
+
+ /* Recurse? */
+ if (c->split) {
+
+ /* Loop over all progeny. Perform pair and self recursion on progenies.*/
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+
+ rec_fof_search_self(props, dim, search_r2, periodic, space_gparts,
+ c->progeny[k]);
+
+ for (int l = k + 1; l < 8; l++)
+ if (c->progeny[l] != NULL)
+ rec_fof_search_pair(props, dim, search_r2, periodic, space_gparts,
+ c->progeny[k], c->progeny[l]);
+ }
+ }
+ }
+ /* Otherwise, compute self-interaction. */
+ else
+ fof_search_self_cell(props, search_r2, space_gparts, c);
+}
+
+/**
+ * @brief Perform the attaching operation using union-find on a given leaf-cell
+ *
+ * @param props The properties fof the FOF scheme.
+ * @param l_x2 The square of the FOF linking length.
+ * @param space_gparts The start of the #gpart array in the #space structure.
+ * @param nr_gparts The number of #gpart in the local #space structure.
+ * @param c The #cell in which to perform FOF.
+ */
+void fof_attach_self_cell(const struct fof_props *props, const double l_x2,
+ const struct gpart *const space_gparts,
+ const size_t nr_gparts, const struct cell *c) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->split) error("Performing the FOF search at a non-leaf level!");
+#endif
+
+ const size_t count = c->grav.count;
+ struct gpart *gparts = (struct gpart *)c->grav.parts;
+
+ /* Make a list of particle offsets into the global gparts array. */
+ size_t *const group_index = props->group_index;
+#ifndef WITH_MPI
+ size_t *const index_offset = group_index + (ptrdiff_t)(gparts - space_gparts);
+#endif
+
+ size_t *const attach_index = props->attach_index;
+ size_t *const attach_offset =
+ attach_index + (ptrdiff_t)(gparts - space_gparts);
+
+ char *const found_attach_index = props->found_attachable_link;
+ char *const found_attach_offset =
+ found_attach_index + (ptrdiff_t)(gparts - space_gparts);
+
+ /* Distances of particles in the global list */
+ float *const offset_dist =
+ props->distance_to_link + (ptrdiff_t)(gparts - space_gparts);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->nodeID != engine_rank)
+ error("Performing self FOF search on foreign cell.");
+#endif
+
+ /* Loop over particles and find which particles belong in the same group. */
+ for (size_t i = 0; i < count; i++) {
+
+ struct gpart *pi = &gparts[i];
+
+ /* Ignore inhibited particles */
+ if (pi->time_bin >= time_bin_inhibited) continue;
+
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(pi)) continue;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pi->ti_drift != ti_current)
+ error("Running FOF on an un-drifted particle!");
+#endif
+
+ const double pix = pi->x[0];
+ const double piy = pi->x[1];
+ const double piz = pi->x[2];
+
+ /* Find the root of pi. */
+#ifdef WITH_MPI
+ const size_t root_i = fof_find_global(
+ i + (ptrdiff_t)(gparts - space_gparts), group_index, nr_gparts);
+#else
+ const size_t root_i = fof_find(index_offset[i], group_index);
+#endif
+
+ /* Get the nature of the linking */
+ const int is_link_i = gpart_is_linkable(pi);
+ const int is_attach_i = gpart_is_attachable(pi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (is_link_i && is_attach_i)
+ error("Particle cannot be both linkable and attachable!");
+#endif
+
+ for (size_t j = i + 1; j < count; j++) {
+
+ struct gpart *pj = &gparts[j];
+
+ /* Ignore inhibited particles */
+ if (pj->time_bin >= time_bin_inhibited) continue;
+
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(pj)) continue;
+
+ /* Get the nature of the linking */
+ const int is_link_j = gpart_is_linkable(pj);
+ const int is_attach_j = gpart_is_attachable(pj);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (is_link_j && is_attach_j)
+ error("Particle cannot be both linkable and attachable!");
+#endif
+
+ /* We only want link<->attach pairs */
+ if (is_attach_i && is_attach_j) continue;
+ if (is_link_i && is_link_j) continue;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pj->ti_drift != ti_current)
+ error("Running FOF on an un-drifted particle!");
+#endif
+
+ /* Find the root of pi. */
+#ifdef WITH_MPI
+ const size_t root_j = fof_find_global(
+ j + (ptrdiff_t)(gparts - space_gparts), group_index, nr_gparts);
+#else
+ const size_t root_j = fof_find(index_offset[j], group_index);
+#endif
+
+ const double pjx = pj->x[0];
+ const double pjy = pj->x[1];
+ const double pjz = pj->x[2];
+
+ /* Compute the pairwise distance */
+ float dx[3], r2 = 0.0f;
+ dx[0] = pix - pjx;
+ dx[1] = piy - pjy;
+ dx[2] = piz - pjz;
+
+ for (int k = 0; k < 3; k++) r2 += dx[k] * dx[k];
+
+ /* Hit or miss? */
+ if (r2 < l_x2) {
+
+ /* Now that we are within the linking length,
+ * decide what to do based on linking types */
+
+ if (is_link_i && is_link_j) {
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Fundamental logic error!");
+#endif
+ } else if (is_link_i && is_attach_j) {
+
+ /* We got a linkable and an attachable.
+ * See whether it is closer and if so re-link.
+ * This is safe to do as the attachables are never roots and
+ * nothing is attached to them */
+ const float dist = sqrtf(r2);
+ if (dist < offset_dist[j]) {
+
+ /* Store the new min dist */
+ offset_dist[j] = dist;
+
+ /* Store the current best root */
+ attach_offset[j] = root_i;
+ found_attach_offset[j] = 1;
+ }
+
+ } else if (is_link_j && is_attach_i) {
+
+ /* We got a linkable and an attachable.
+ * See whether it is closer and if so re-link.
+ * This is safe to do as the attachables are never roots and
+ * nothing is attached to them */
+ const float dist = sqrtf(r2);
+ if (dist < offset_dist[i]) {
+
+ /* Store the new min dist */
+ offset_dist[i] = dist;
+
+ /* Store the current best root */
+ attach_offset[i] = root_j;
+ found_attach_offset[i] = 1;
+ }
+
+ } else {
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Fundamental logic error!");
+#endif
+ }
+ }
+ }
+ }
+}
+
+/**
+ * @brief Perform the attaching operation using union-find between two cells
+ *
+ * @param props The properties fof the FOF scheme.
+ * @param dim The dimension of the simulation volume.
+ * @param l_x2 The square of the FOF linking length.
+ * @param periodic Are we using periodic BCs?
+ * @param space_gparts The start of the #gpart array in the #space structure.
+ * @param nr_gparts The number of #gpart in the local #space structure.
+ * @param ci The first #cell in which to perform FOF.
+ * @param cj The second #cell in which to perform FOF.
+ * @param ci_local Is the #cell ci on the local MPI rank?
+ * @param cj_local Is the #cell cj on the local MPI rank?
+ */
+void fof_attach_pair_cells(const struct fof_props *props, const double dim[3],
+ const double l_x2, const int periodic,
+ const struct gpart *const space_gparts,
+ const size_t nr_gparts,
+ const struct cell *restrict ci,
+ const struct cell *restrict cj, const int ci_local,
+ const int cj_local) {
+
+ const size_t count_i = ci->grav.count;
+ const size_t count_j = cj->grav.count;
+ struct gpart *gparts_i = (struct gpart *)ci->grav.parts;
+ struct gpart *gparts_j = (struct gpart *)cj->grav.parts;
+
+ /* Index of particles in the global group list */
+ size_t *const group_index = props->group_index;
+
+ /* Make a list of particle offsets into the global gparts array. */
+ size_t *const index_offset_i =
+ group_index + (ptrdiff_t)(gparts_i - space_gparts);
+ size_t *const index_offset_j =
+ group_index + (ptrdiff_t)(gparts_j - space_gparts);
+
+ size_t *const attach_offset_i =
+ props->attach_index + (ptrdiff_t)(gparts_i - space_gparts);
+ size_t *const attach_offset_j =
+ props->attach_index + (ptrdiff_t)(gparts_j - space_gparts);
+
+ char *const found_attach_offset_i =
+ props->found_attachable_link + (ptrdiff_t)(gparts_i - space_gparts);
+ char *const found_attach_offset_j =
+ props->found_attachable_link + (ptrdiff_t)(gparts_j - space_gparts);
+
+ /* Distances of particles in the global list */
+ float *const offset_dist_i =
+ props->distance_to_link + (ptrdiff_t)(gparts_i - space_gparts);
+ float *const offset_dist_j =
+ props->distance_to_link + (ptrdiff_t)(gparts_j - space_gparts);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (index_offset_j > index_offset_i &&
+ (index_offset_j < index_offset_i + count_i))
+ error("Overlapping cells");
+ if (index_offset_i > index_offset_j &&
+ (index_offset_i < index_offset_j + count_j))
+ error("Overlapping cells");
+#endif
+
+ /* Account for boundary conditions.*/
+ double shift[3] = {0.0, 0.0, 0.0};
+
+ /* Get the relative distance between the pairs, wrapping. */
+ double diff[3];
+ for (int k = 0; k < 3; k++) {
+ diff[k] = cj->loc[k] - ci->loc[k];
+ if (periodic && diff[k] < -dim[k] * 0.5)
+ shift[k] = dim[k];
+ else if (periodic && diff[k] > dim[k] * 0.5)
+ shift[k] = -dim[k];
+ else
+ shift[k] = 0.0;
+ diff[k] += shift[k];
+ }
+
+ /* Loop over particles and find which particles belong in the same group. */
+ for (size_t i = 0; i < count_i; i++) {
+
+ struct gpart *restrict pi = &gparts_i[i];
+
+ /* Ignore inhibited particles */
+ if (pi->time_bin >= time_bin_inhibited) continue;
+
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(pi)) continue;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pi->ti_drift != ti_current)
+ error("Running FOF on an un-drifted particle!");
+#endif
+
+ const double pix = pi->x[0] - shift[0];
+ const double piy = pi->x[1] - shift[1];
+ const double piz = pi->x[2] - shift[2];
+
+ /* Find the root of pi. */
+#ifdef WITH_MPI
+ size_t root_i;
+ if (ci_local) {
+ root_i = fof_find_global(index_offset_i[i] - node_offset, group_index,
+ nr_gparts);
+ } else {
+ root_i = pi->fof_data.group_id;
+ }
+#else
+ const size_t root_i = fof_find(index_offset_i[i], group_index);
+#endif
+
+ /* Get the nature of the linking */
+ const int is_link_i = gpart_is_linkable(pi);
+ const int is_attach_i = gpart_is_attachable(pi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (is_link_i && is_attach_i)
+ error("Particle cannot be both linkable and attachable!");
+#endif
+
+ for (size_t j = 0; j < count_j; j++) {
+
+ struct gpart *restrict pj = &gparts_j[j];
+
+ /* Ignore inhibited particles */
+ if (pj->time_bin >= time_bin_inhibited) continue;
+
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(pj)) continue;
+
+ /* Get the nature of the linking */
+ const int is_link_j = gpart_is_linkable(pj);
+ const int is_attach_j = gpart_is_attachable(pj);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (is_link_j && is_attach_j)
+ error("Particle cannot be both linkable and attachable!");
+#endif
+
+ /* We only want link<->attach pairs */
+ if (is_attach_i && is_attach_j) continue;
+ if (is_link_i && is_link_j) continue;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pj->ti_drift != ti_current)
+ error("Running FOF on an un-drifted particle!");
+#endif
+
+ /* Find the root of pj. */
+#ifdef WITH_MPI
+ size_t root_j;
+ if (cj_local) {
+ root_j = fof_find_global(index_offset_j[j] - node_offset, group_index,
+ nr_gparts);
+ } else {
+ root_j = pj->fof_data.group_id;
+ }
+#else
+ const size_t root_j = fof_find(index_offset_j[j], group_index);
+#endif
+
+ const double pjx = pj->x[0];
+ const double pjy = pj->x[1];
+ const double pjz = pj->x[2];
+
+ /* Compute pairwise distance (periodic BCs were accounted
+ for by the shift vector) */
+ float dx[3], r2 = 0.0f;
+ dx[0] = pix - pjx;
+ dx[1] = piy - pjy;
+ dx[2] = piz - pjz;
+
+ for (int k = 0; k < 3; k++) r2 += dx[k] * dx[k];
+
+ /* Hit or miss? */
+ if (r2 < l_x2) {
+
+ /* Now that we are within the linking length,
+ * decide what to do based on linking types */
+
+ if (is_link_i && is_link_j) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Fundamental logic error!");
+#endif
+
+ } else if (is_link_i && is_attach_j) {
+
+ /* We got a linkable and an attachable.
+ * See whether it is closer and if so re-link.
+ * This is safe to do as the attachables are never roots and
+ * nothing is attached to them */
+ if (cj_local) {
+
+ const float dist = sqrtf(r2);
+
+ if (dist < offset_dist_j[j]) {
+
+ /* Store the new min dist */
+ offset_dist_j[j] = dist;
+
+ /* Store the current best root */
+ attach_offset_j[j] = root_i;
+ found_attach_offset_j[j] = 1;
+ }
+ }
+
+ } else if (is_link_j && is_attach_i) {
+
+ /* We got a linkable and an attachable.
+ * See whether it is closer and if so re-link.
+ * This is safe to do as the attachables are never roots and
+ * nothing is attached to them */
+ if (ci_local) {
+
+ const float dist = sqrtf(r2);
+
+ if (dist < offset_dist_i[i]) {
+
+ /* Store the new min dist */
+ offset_dist_i[i] = dist;
+
+ /* Store the current best root */
+ attach_offset_i[i] = root_j;
+ found_attach_offset_i[i] = 1;
+ }
+ }
+
+ } else {
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Fundamental logic error!");
+#endif
+ }
+ }
+ }
+ }
+}
+
+/**
+ * @brief Recursively perform a union-find attaching between two cells.
+ *
+ * If cells are more distant than the linking length, we abort early.
+ *
+ * @param props The properties fof the FOF scheme.
+ * @param dim The dimension of the space.
+ * @param attach_r2 the square of the FOF linking length.
+ * @param periodic Are we using periodic BCs?
+ * @param space_gparts The start of the #gpart array in the #space structure.
+ * @param nr_gparts The number of #gpart in the local #space structure.
+ * @param ci The first #cell in which to perform FOF.
+ * @param cj The second #cell in which to perform FOF.
+ * @param ci_local Is the #cell ci on the local MPI rank?
+ * @param cj_local Is the #cell cj on the local MPI rank?
+ */
+void rec_fof_attach_pair(const struct fof_props *props, const double dim[3],
+ const double attach_r2, const int periodic,
+ const struct gpart *const space_gparts,
+ const size_t nr_gparts, struct cell *restrict ci,
+ struct cell *restrict cj, const int ci_local,
+ const int cj_local) {
/* Find the shortest distance between cells, remembering to account for
* boundary conditions. */
const double r2 = cell_min_dist(ci, cj, dim);
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci == cj) error("Pair FOF called on same cell!!!");
+#endif
+
/* Return if cells are out of range of each other. */
- if (r2 > search_r2) return;
+ if (r2 > attach_r2) return;
/* Recurse on both cells if they are both split. */
if (ci->split && cj->split) {
@@ -1244,52 +1994,46 @@ void rec_fof_search_pair_foreign(
for (int l = 0; l < 8; l++)
if (cj->progeny[l] != NULL)
- rec_fof_search_pair_foreign(props, dim, search_r2, periodic,
- space_gparts, nr_gparts, ci->progeny[k],
- cj->progeny[l], link_count, group_links,
- group_links_size);
+ rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts,
+ nr_gparts, ci->progeny[k], cj->progeny[l],
+ ci_local, cj_local);
}
}
} else if (ci->split) {
-
for (int k = 0; k < 8; k++) {
if (ci->progeny[k] != NULL)
- rec_fof_search_pair_foreign(props, dim, search_r2, periodic,
- space_gparts, nr_gparts, ci->progeny[k], cj,
- link_count, group_links, group_links_size);
+ rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts,
+ nr_gparts, ci->progeny[k], cj, ci_local, cj_local);
}
} else if (cj->split) {
for (int k = 0; k < 8; k++) {
if (cj->progeny[k] != NULL)
- rec_fof_search_pair_foreign(props, dim, search_r2, periodic,
- space_gparts, nr_gparts, ci, cj->progeny[k],
- link_count, group_links, group_links_size);
+ rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts,
+ nr_gparts, ci, cj->progeny[k], ci_local, cj_local);
}
} else {
- /* Perform FOF search between pairs of cells that are within the linking
+ /* Perform FOF attach between pairs of cells that are within the linking
* length and not the same cell. */
- fof_search_pair_cells_foreign(props, dim, search_r2, periodic, space_gparts,
- nr_gparts, ci, cj, link_count, group_links,
- group_links_size);
+ fof_attach_pair_cells(props, dim, attach_r2, periodic, space_gparts,
+ nr_gparts, ci, cj, ci_local, cj_local);
}
}
-#endif /* WITH_MPI */
-
/**
- * @brief Recursively perform a union-find FOF on a cell.
+ * @brief Recursively perform a the attaching operation on a cell.
*
* @param props The properties fof the FOF scheme.
* @param dim The dimension of the space.
- * @param space_gparts The start of the #gpart array in the #space structure.
- * @param search_r2 the square of the FOF linking length.
+ * @param attach_r2 the square of the FOF linking length.
* @param periodic Are we using periodic BCs?
+ * @param space_gparts The start of the #gpart array in the #space structure.
+ * @param nr_gparts The number of #gpart in the local #space structure.
* @param c The #cell in which to perform FOF.
*/
-void rec_fof_search_self(const struct fof_props *props, const double dim[3],
- const double search_r2, const int periodic,
+void rec_fof_attach_self(const struct fof_props *props, const double dim[3],
+ const double attach_r2, const int periodic,
const struct gpart *const space_gparts,
- struct cell *c) {
+ const size_t nr_gparts, struct cell *c) {
/* Recurse? */
if (c->split) {
@@ -1298,19 +2042,22 @@ void rec_fof_search_self(const struct fof_props *props, const double dim[3],
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
- rec_fof_search_self(props, dim, search_r2, periodic, space_gparts,
- c->progeny[k]);
+ rec_fof_attach_self(props, dim, attach_r2, periodic, space_gparts,
+ nr_gparts, c->progeny[k]);
for (int l = k + 1; l < 8; l++)
if (c->progeny[l] != NULL)
- rec_fof_search_pair(props, dim, search_r2, periodic, space_gparts,
- c->progeny[k], c->progeny[l]);
+ rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts,
+ nr_gparts, c->progeny[k], c->progeny[l],
+ /*ci_local=*/1,
+ /*cj_local=*/1);
}
}
+ } else {
+
+ /* Otherwise, compute self-interaction. */
+ fof_attach_self_cell(props, attach_r2, space_gparts, nr_gparts, c);
}
- /* Otherwise, compute self-interaction. */
- else
- fof_search_self_cell(props, search_r2, space_gparts, c);
}
/* Mapper function to atomically update the group size array. */
@@ -1336,22 +2083,20 @@ void fof_calc_group_size_mapper(void *map_data, int num_elements,
/* Retrieve mapped data. */
struct space *s = (struct space *)extra_data;
struct gpart *gparts = (struct gpart *)map_data;
- size_t *group_index = s->e->fof_properties->group_index;
- size_t *group_size = s->e->fof_properties->group_size;
+ size_t *restrict group_index = s->e->fof_properties->group_index;
+ size_t *restrict group_size = s->e->fof_properties->group_size;
/* Offset into gparts array. */
- ptrdiff_t gparts_offset = (ptrdiff_t)(gparts - s->gparts);
+ const ptrdiff_t gparts_offset = (ptrdiff_t)(gparts - s->gparts);
size_t *const group_index_offset = group_index + gparts_offset;
/* Create hash table. */
hashmap_t map;
hashmap_init(&map);
- /* Loop over particles and find which cells are in range of each other to
- * perform the FOF search. */
for (int ind = 0; ind < num_elements; ind++) {
- hashmap_key_t root =
+ const hashmap_key_t root =
(hashmap_key_t)fof_find(group_index_offset[ind], group_index);
const size_t gpart_index = gparts_offset + ind;
@@ -1375,9 +2120,9 @@ void fof_calc_group_size_mapper(void *map_data, int num_elements,
}
/* Mapper function to atomically update the group mass array. */
-static INLINE void fof_update_group_mass_mapper(hashmap_key_t key,
- hashmap_value_t *value,
- void *data) {
+static INLINE void fof_update_group_mass_iterator(hashmap_key_t key,
+ hashmap_value_t *value,
+ void *data) {
double *group_mass = (double *)data;
@@ -1385,6 +2130,16 @@ static INLINE void fof_update_group_mass_mapper(hashmap_key_t key,
atomic_add_d(&group_mass[key], value->value_dbl);
}
+/* Mapper function to atomically update the group size array. */
+static INLINE void fof_update_group_size_iterator(hashmap_key_t key,
+ hashmap_value_t *value,
+ void *data) {
+ long long *group_size = (long long *)data;
+
+ /* Use key to index into group mass array. */
+ atomic_add(&group_size[key], value->value_st);
+}
+
/**
* @brief Mapper function to calculate the group masses.
*
@@ -1399,6 +2154,7 @@ void fof_calc_group_mass_mapper(void *map_data, int num_elements,
struct space *s = (struct space *)extra_data;
struct gpart *gparts = (struct gpart *)map_data;
double *group_mass = s->e->fof_properties->group_mass;
+ long long *group_size = s->e->fof_properties->final_group_size;
const size_t group_id_default = s->e->fof_properties->group_id_default;
const size_t group_id_offset = s->e->fof_properties->group_id_offset;
@@ -1417,16 +2173,19 @@ void fof_calc_group_mass_mapper(void *map_data, int num_elements,
hashmap_value_t *data = hashmap_get(&map, index);
/* Update group mass */
- if (data != NULL)
+ if (data != NULL) {
(*data).value_dbl += gparts[ind].mass;
- else
+ (*data).value_st += 1;
+ } else
error("Couldn't find key (%zu) or create new one.", index);
}
}
/* Update the group mass array. */
if (map.size > 0)
- hashmap_iterate(&map, fof_update_group_mass_mapper, group_mass);
+ hashmap_iterate(&map, fof_update_group_mass_iterator, group_mass);
+ if (map.size > 0)
+ hashmap_iterate(&map, fof_update_group_size_iterator, group_size);
hashmap_free(&map);
}
@@ -1444,6 +2203,7 @@ void fof_unpack_group_mass_mapper(hashmap_key_t key, hashmap_value_t *value,
/* Store elements from hash table in array. */
mass_send[*nsend].global_root = key;
mass_send[*nsend].group_mass = value->value_dbl;
+ mass_send[*nsend].final_group_size = value->value_ll;
mass_send[*nsend].first_position[0] = value->value_array2_dbl[0];
mass_send[*nsend].first_position[1] = value->value_array2_dbl[1];
mass_send[*nsend].first_position[2] = value->value_array2_dbl[2];
@@ -1488,6 +2248,7 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
float *max_part_density = props->max_part_density;
double *centre_of_mass = props->group_centre_of_mass;
double *first_position = props->group_first_position;
+ long long *final_group_size = props->final_group_size;
/* Start the hash map */
hashmap_t map;
@@ -1500,8 +2261,8 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
/* Ignore inhibited particles */
if (gparts[i].time_bin >= time_bin_inhibited) continue;
- /* Ignore neutrinos */
- if (gparts[i].type == swift_type_neutrino) continue;
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(&gparts[i])) continue;
/* Check if the particle is in a group above the threshold. */
if (gparts[i].fof_data.group_id != group_id_default) {
@@ -1518,6 +2279,9 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
/* Update group mass */
group_mass[index] += gparts[i].mass;
+ /* Update group size */
+ final_group_size[index]++;
+
} else {
/* The root is *not* local */
@@ -1534,6 +2298,9 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
/* Add mass fragments of groups */
data->value_dbl += mass;
+ /* Increase fragment size */
+ data->value_ll++;
+
/* Record the first particle of this fragment that we encounter so we
* we can use it as reference frame for the centre of mass calculation
*/
@@ -1648,6 +2415,7 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
const size_t index =
gparts[local_root_index].fof_data.group_id - local_group_offset;
group_mass[index] += fof_mass_recv[i].group_mass;
+ final_group_size[index] += fof_mass_recv[i].final_group_size;
}
/* Loop over particles, densest particle in each *local* group.
@@ -1657,8 +2425,8 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
/* Ignore inhibited particles */
if (gparts[i].time_bin >= time_bin_inhibited) continue;
- /* Ignore neutrinos */
- if (gparts[i].type == swift_type_neutrino) continue;
+ /* Check whether we ignore this particle type altogether */
+ if (current_fof_ignore_type & (1 << (gparts[i].type + 1))) continue;
/* Only check groups above the minimum mass threshold. */
if (gparts[i].fof_data.group_id != group_id_default) {
@@ -1879,8 +2647,8 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
/* Ignore inhibited particles */
if (gparts[i].time_bin >= time_bin_inhibited) continue;
- /* Ignore neutrinos */
- if (gparts[i].type == swift_type_neutrino) continue;
+ /* Check whether we ignore this particle type altogether */
+ if (gpart_is_ignorable(&gparts[i])) continue;
const size_t index = gparts[i].fof_data.group_id - group_id_offset;
@@ -1982,8 +2750,8 @@ void fof_find_foreign_links_mapper(void *map_data, int num_elements,
for (int ind = 0; ind < num_elements; ind++) {
/* Get the local and foreign cells to recurse on. */
- struct cell *restrict local_cell = cell_pairs[ind].local;
- struct cell *restrict foreign_cell = cell_pairs[ind].foreign;
+ const struct cell *restrict local_cell = cell_pairs[ind].local;
+ const struct cell *restrict foreign_cell = cell_pairs[ind].foreign;
rec_fof_search_pair_foreign(props, dim, search_r2, periodic, gparts,
nr_gparts, local_cell, foreign_cell,
@@ -1996,8 +2764,8 @@ void fof_find_foreign_links_mapper(void *map_data, int num_elements,
if (lock_lock(&s->lock) == 0) {
/* Get pointers to global arrays. */
- int *group_links_size = &props->group_links_size;
- int *group_link_count = &props->group_link_count;
+ int *restrict group_links_size = &props->group_links_size;
+ int *restrict group_link_count = &props->group_link_count;
struct fof_mpi **group_links = &props->group_links;
/* If the global group_links array is not big enough re-allocate it. */
@@ -2244,7 +3012,7 @@ void fof_dump_group_data(const struct fof_props *props, const int my_rank,
FILE *file = NULL;
struct part *parts = s->parts;
- size_t *group_size = props->group_size;
+ long long *final_group_size = props->final_group_size;
size_t *group_index = props->group_index;
double *group_mass = props->group_mass;
double *group_centre_of_mass = props->group_centre_of_mass;
@@ -2287,8 +3055,8 @@ void fof_dump_group_data(const struct fof_props *props, const int my_rank,
const long long part_id = props->max_part_density_index[i] >= 0
? parts[max_part_density_index[i]].id
: -1;
- fprintf(file, " %8zu %12zu %12e %12e %12e %12e %12e %24lld %24lld\n",
- group_index[i], group_size[i], group_mass[i],
+ fprintf(file, " %8zu %12lld %12e %12e %12e %12e %12e %24lld %24lld\n",
+ group_index[i], final_group_size[i], group_mass[i],
group_centre_of_mass[i * 3 + 0],
group_centre_of_mass[i * 3 + 1],
group_centre_of_mass[i * 3 + 2], max_part_density[i],
@@ -2313,6 +3081,7 @@ void fof_dump_group_data(const struct fof_props *props, const int my_rank,
struct mapper_data {
size_t *group_index;
size_t *group_size;
+ float *distance_to_link;
size_t nr_gparts;
struct gpart *space_gparts;
};
@@ -2349,9 +3118,19 @@ void fof_set_outgoing_root_mapper(void *map_data, int num_elements,
/* Set each particle's root and group properties found in the local FOF.*/
for (int k = 0; k < local_cell->grav.count; k++) {
+
+ /* TODO: Can we skip ignorable particles here?
+ * Likely makes no difference */
+
+ /* Recall we did alter the group_index with a global_offset.
+ * We need to remove that here as we want the *local* root */
const size_t root =
fof_find_global(offset[k] - node_offset, group_index, nr_gparts);
+ /* TODO: Could we call fof_find() here instead?
+ * Likely yes but we don't want path compression at this stage.
+ * So, probably not */
+
gparts[k].fof_data.group_id = root;
gparts[k].fof_data.group_size = group_size[root - node_offset];
}
@@ -2372,31 +3151,33 @@ void fof_set_outgoing_root_mapper(void *map_data, int num_elements,
void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
#ifdef WITH_MPI
-
struct engine *e = s->e;
const int verbose = e->verbose;
+
+ /* Abort if only one node */
+ if (e->nr_nodes == 1) return;
+
size_t *restrict group_index = props->group_index;
size_t *restrict group_size = props->group_size;
const size_t nr_gparts = s->nr_gparts;
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const double search_r2 = props->l_x2;
+ const ticks tic_total = getticks();
ticks tic = getticks();
/* Make group IDs globally unique. */
for (size_t i = 0; i < nr_gparts; i++) group_index[i] += node_offset;
struct cell_pair_indices *cell_pairs = NULL;
- int group_link_count = 0;
int cell_pair_count = 0;
props->group_links_size = fof_props_default_group_link_size;
- props->group_link_count = 0;
int num_cells_out = 0;
int num_cells_in = 0;
- /* Find the maximum no. of cell pairs. */
+ /* Find the maximum no. of cell pairs that can communicate. */
for (int i = 0; i < e->nr_proxies; i++) {
for (int j = 0; j < e->proxies[i].nr_cells_out; j++) {
@@ -2422,6 +3203,7 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
const int cell_pair_size = num_cells_in * num_cells_out;
+ /* Allocate memory for all the possible cell links */
if (swift_memalign("fof_group_links", (void **)&props->group_links,
SWIFT_STRUCT_ALIGNMENT,
props->group_links_size * sizeof(struct fof_mpi)) != 0)
@@ -2464,12 +3246,13 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
/* Skip empty cells. */
if (foreign_cell->grav.count == 0) continue;
- /* Check if local cell has already been added to the local list of
- * cells. */
+ /* Add candidates in range to the list of pairs of cells to treat. */
const double r2 = cell_min_dist(local_cell, foreign_cell, dim);
if (r2 < search_r2) {
cell_pairs[cell_pair_count].local = local_cell;
- cell_pairs[cell_pair_count++].foreign = foreign_cell;
+ cell_pairs[cell_pair_count].foreign = foreign_cell;
+
+ cell_pair_count++;
}
}
}
@@ -2477,12 +3260,10 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
}
if (verbose)
- message(
- "Finding local/foreign cell pairs"
- "took: %.3f %s.",
- clocks_from_ticks(getticks() - tic_pairs), clocks_getunit());
+ message("Finding local/foreign cell pairs took: %.3f %s.",
+ clocks_from_ticks(getticks() - tic_pairs), clocks_getunit());
- ticks tic_set_roots = getticks();
+ const ticks tic_set_roots = getticks();
/* Set the root of outgoing particles. */
@@ -2502,7 +3283,7 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
}
}
- /* Now set the roots */
+ /* Now set the *local* roots of all the gparts we are sending */
struct mapper_data data;
data.group_index = group_index;
data.group_size = group_size;
@@ -2513,10 +3294,8 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
threadpool_auto_chunk_size, &data);
if (verbose)
- message(
- "Initialising particle roots "
- "took: %.3f %s.",
- clocks_from_ticks(getticks() - tic_set_roots), clocks_getunit());
+ message("Initialising particle roots took: %.3f %s.",
+ clocks_from_ticks(getticks() - tic_set_roots), clocks_getunit());
free(local_cells);
@@ -2526,14 +3305,12 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
"took: %.3f %s.",
clocks_from_ticks(getticks() - tic), clocks_getunit());
- /* Allocate buffers to receive the gpart fof information */
- engine_allocate_foreign_particles(e, /*fof=*/1);
-
/* Activate the tasks exchanging all the required gparts */
engine_activate_gpart_comms(e);
ticks local_fof_tic = getticks();
+ /* Wait for all the communication tasks to be ready */
MPI_Barrier(MPI_COMM_WORLD);
if (verbose)
@@ -2549,19 +3326,22 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
message("MPI send/recv comms took: %.3f %s.",
clocks_from_ticks(getticks() - tic), clocks_getunit());
+ /* We have now recevied the foreign particles. Each particle received
+ * carries information about its own *foreign* (to us) root and the
+ * size of the group fragment it belongs too its original foreign rank. */
+
tic = getticks();
+ props->group_link_count = 0;
+
/* Perform search of group links between local and foreign cells with the
* threadpool. */
threadpool_map(&s->e->threadpool, fof_find_foreign_links_mapper, cell_pairs,
cell_pair_count, sizeof(struct cell_pair_indices), 1,
(struct space *)s);
- group_link_count = props->group_link_count;
-
- /* Clean up memory. */
+ /* Clean up memory used by foreign particles. */
swift_free("fof_cell_pairs", cell_pairs);
- space_free_foreign_parts(e->s, /*clear pointers=*/1);
if (verbose)
message("Searching for foreign links took: %.3f %s.",
@@ -2569,11 +3349,7 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
tic = getticks();
- struct fof_mpi *global_group_links = NULL;
- int *displ = NULL, *group_link_counts = NULL;
- int global_group_link_count = 0;
-
- ticks comms_tic = getticks();
+ const ticks comms_tic = getticks();
MPI_Barrier(MPI_COMM_WORLD);
@@ -2581,14 +3357,168 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
message("Imbalance took: %.3f %s.",
clocks_from_ticks(getticks() - comms_tic), clocks_getunit());
- comms_tic = getticks();
+ if (verbose)
+ message("fof_search_foreign_cells() took (FOF SCALING): %.3f %s.",
+ clocks_from_ticks(getticks() - tic_total), clocks_getunit());
+
+#endif /* WITH_MPI */
+}
+
+/**
+ * @brief Run all the tasks attaching the attachables to their
+ * nearest linkable particle.
+ *
+ * @param props The properties fof the FOF scheme.
+ * @param s The #space we work with.
+ */
+void fof_link_attachable_particles(struct fof_props *props,
+ const struct space *s) {
+
+ /* Is there anything to attach? */
+ if (!current_fof_attach_type) return;
+
+ const ticks tic_total = getticks();
+
+ /* Activate the tasks attaching attachable particles to the linkable ones */
+ engine_activate_fof_attach_tasks(s->e);
+
+ /* Perform FOF tasks for attachable particles. */
+ engine_launch(s->e, "fof");
+
+ if (s->e->verbose)
+ message("fof_link_attachable_particles() took (FOF SCALING): %.3f %s.",
+ clocks_from_ticks(getticks() - tic_total), clocks_getunit());
+}
+
+/**
+ * @brief Process all the attachable-linkable connections to add the
+ * attachables to the groups they belong to.
+ *
+ * @param props The properties fof the FOF scheme.
+ * @param s The #space we work with.
+ */
+void fof_finalise_attachables(struct fof_props *props, const struct space *s) {
+
+ /* Is there anything to attach? */
+ if (!current_fof_attach_type) return;
+
+ const ticks tic_total = getticks();
+
+ const size_t nr_gparts = s->nr_gparts;
+
+ char *found_attachable_link = props->found_attachable_link;
+ size_t *restrict group_index = props->group_index;
+ size_t *restrict attach_index = props->attach_index;
+ size_t *restrict group_size = props->group_size;
+
+#ifdef WITH_MPI
+
+ /* Get pointers to global arrays. */
+ int *restrict group_links_size = &props->group_links_size;
+ int *restrict group_link_count = &props->group_link_count;
+ struct fof_mpi **group_links = &props->group_links;
+
+ /* Loop over all the attachables and added them to the group they belong to */
+ for (size_t i = 0; i < nr_gparts; ++i) {
+
+ const struct gpart *gp = &s->gparts[i];
+
+ if (gpart_is_attachable(gp) && found_attachable_link[i]) {
+
+ /* Update its root */
+ const size_t root_j = attach_index[i];
+ const size_t root_i =
+ fof_find_global(group_index[i] - node_offset, group_index, nr_gparts);
+
+ /* Update the size of the group the particle belongs to */
+ if (is_local(root_j, nr_gparts)) {
+
+ const size_t local_root = root_j - node_offset;
+
+ group_index[i] = local_root + node_offset;
+ group_size[local_root]++;
+
+ } else {
+
+ add_foreign_link_to_list(group_link_count, group_links_size,
+ group_links, group_links, root_i, root_j,
+ /*size_i=*/1,
+ /*size_j=*/2);
+ }
+ }
+ }
+
+#else /* not WITH_MPI */
+
+ /* Loop over all the attachables and added them to the group they belong to */
+ for (size_t i = 0; i < nr_gparts; ++i) {
+
+ const struct gpart *gp = &s->gparts[i];
+
+ if (gpart_is_attachable(gp) && found_attachable_link[i]) {
+
+ const size_t root = attach_index[i];
+
+ /* Update its root */
+ group_index[i] = root;
+
+ /* Update the size of the group the particle belongs to */
+ group_size[root]++;
+ }
+ }
+
+#endif /* WITH_MPI */
+
+ if (s->e->verbose)
+ message("fof_finalise_attachables() took (FOF SCALING): %.3f %s.",
+ clocks_from_ticks(getticks() - tic_total), clocks_getunit());
+}
+
+/**
+ * @brief Process all the group fragments spanning more than
+ * one rank to link them.
+ *
+ * This is the final global union-find pass which concludes
+ * the MPI-FOF-algorithm.
+ *
+ * @param props The properties fof the FOF scheme.
+ * @param s The #space we work with.
+ */
+void fof_link_foreign_fragments(struct fof_props *props,
+ const struct space *s) {
+
+#ifdef WITH_MPI
+
+ struct engine *e = s->e;
+ const int verbose = e->verbose;
+
+ /* Abort if only one node */
+ if (e->nr_nodes == 1) return;
+
+ const size_t nr_gparts = s->nr_gparts;
+ size_t *restrict group_index = props->group_index;
+ size_t *restrict group_size = props->group_size;
+
+ const ticks tic_total = getticks();
+ ticks tic = getticks();
+ const ticks comms_tic = getticks();
+
+ if (verbose)
+ message(
+ "Searching %zu gravity particles for cross-node links with l_x: %lf",
+ nr_gparts, sqrt(props->l_x2));
+
+ /* Local copy of the variable set in the mapper */
+ const int group_link_count = props->group_link_count;
/* Sum the total number of links across MPI domains over each MPI rank. */
+ int global_group_link_count = 0;
MPI_Allreduce(&group_link_count, &global_group_link_count, 1, MPI_INT,
MPI_SUM, MPI_COMM_WORLD);
- /* Unique set of links is half of all group links as each link is found twice
- * by opposing MPI ranks. */
+ struct fof_mpi *global_group_links = NULL;
+ int *displ = NULL, *group_link_counts = NULL;
+
if (swift_memalign("fof_global_group_links", (void **)&global_group_links,
SWIFT_STRUCT_ALIGNMENT,
global_group_link_count * sizeof(struct fof_mpi)) != 0)
@@ -2613,8 +3543,10 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
/* Set the displacements into the global link list using the link counts from
* each rank */
displ[0] = 0;
- for (int i = 1; i < e->nr_nodes; i++)
+ for (int i = 1; i < e->nr_nodes; i++) {
displ[i] = displ[i - 1] + group_link_counts[i - 1];
+ if (displ[i] < 0) error("Number of group links overflowing!");
+ }
/* Gather the global link list on all ranks. */
MPI_Allgatherv(props->group_links, group_link_count, fof_mpi_type,
@@ -2638,11 +3570,11 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
tic = getticks();
/* Transform the group IDs to a local list going from 0-group_count so a
- * union-find can be performed. */
+ * union-find can be performed.
+ * Each member of a link is stored separately --> Need 2x as many entries */
size_t *global_group_index = NULL, *global_group_id = NULL,
*global_group_size = NULL;
const int global_group_list_size = 2 * global_group_link_count;
- int group_count = 0;
if (swift_memalign("fof_global_group_index", (void **)&global_group_index,
SWIFT_STRUCT_ALIGNMENT,
@@ -2672,18 +3604,21 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
hashmap_init(&map);
/* Store each group ID and its properties. */
- for (int i = 0; i < global_group_link_count; i++) {
+ int group_count = 0;
+ for (int k = 0; k < global_group_link_count; k++) {
- size_t group_i = global_group_links[i].group_i;
- size_t group_j = global_group_links[i].group_j;
+ const size_t group_i = global_group_links[k].group_i;
+ const size_t group_j = global_group_links[k].group_j;
- global_group_size[group_count] += global_group_links[i].group_i_size;
+ global_group_size[group_count] += global_group_links[k].group_i_size;
global_group_id[group_count] = group_i;
- hashmap_add_group(group_i, group_count++, &map);
+ hashmap_add_group(group_i, group_count, &map);
+ group_count++;
- global_group_size[group_count] += global_group_links[i].group_j_size;
+ global_group_size[group_count] += global_group_links[k].group_j_size;
global_group_id[group_count] = group_j;
- hashmap_add_group(group_j, group_count++, &map);
+ hashmap_add_group(group_j, group_count, &map);
+ group_count++;
}
if (verbose)
@@ -2693,8 +3628,8 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
tic = getticks();
/* Create a global_group_index list of groups across MPI domains so that you
- * can perform a union-find locally on each node. */
- /* The value of which is an offset into global_group_id, which is the actual
+ * can perform a union-find locally on each node.
+ * The value of which is an offset into global_group_id, which is the actual
* root. */
for (int i = 0; i < group_count; i++) global_group_index[i] = i;
@@ -2708,30 +3643,30 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
"Error while allocating memory for the displacement in memory for the "
"global group link list");
- for (int i = 0; i < group_count; i++)
- orig_global_group_size[i] = global_group_size[i];
+ memcpy(orig_global_group_size, global_group_size,
+ group_count * sizeof(size_t));
/* Perform a union-find on the group links. */
- for (int i = 0; i < global_group_link_count; i++) {
+ for (int k = 0; k < global_group_link_count; k++) {
/* Use the hash table to find the group offsets in the index array. */
- size_t find_i =
- hashmap_find_group_offset(global_group_links[i].group_i, &map);
- size_t find_j =
- hashmap_find_group_offset(global_group_links[i].group_j, &map);
+ const size_t find_i =
+ hashmap_find_group_offset(global_group_links[k].group_i, &map);
+ const size_t find_j =
+ hashmap_find_group_offset(global_group_links[k].group_j, &map);
/* Use the offset to find the group's root. */
- size_t root_i = fof_find(find_i, global_group_index);
- size_t root_j = fof_find(find_j, global_group_index);
+ const size_t root_i = fof_find(find_i, global_group_index);
+ const size_t root_j = fof_find(find_j, global_group_index);
- size_t group_i = global_group_id[root_i];
- size_t group_j = global_group_id[root_j];
+ const size_t group_i = global_group_id[root_i];
+ const size_t group_j = global_group_id[root_j];
if (group_i == group_j) continue;
/* Update roots accordingly. */
- size_t size_i = global_group_size[root_i];
- size_t size_j = global_group_size[root_j];
+ const size_t size_i = global_group_size[root_i];
+ const size_t size_j = global_group_size[root_j];
#ifdef UNION_BY_SIZE_OVER_MPI
if (size_i < size_j) {
global_group_index[root_i] = root_j;
@@ -2762,9 +3697,9 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
/* Update each group locally with new root information. */
for (int i = 0; i < group_count; i++) {
- size_t group_id = global_group_id[i];
- size_t offset = fof_find(global_group_index[i], global_group_index);
- size_t new_root = global_group_id[offset];
+ const size_t group_id = global_group_id[i];
+ const size_t offset = fof_find(global_group_index[i], global_group_index);
+ const size_t new_root = global_group_id[offset];
/* If the group is local update its root and size. */
if (is_local(group_id, nr_gparts) && new_root != group_id) {
@@ -2792,67 +3727,39 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
swift_free("fof_global_group_id", global_group_id);
swift_free("fof_orig_global_group_size", orig_global_group_size);
+ if (verbose) {
+ message("link_foreign_fragmens() took (FOF SCALING): %.3f %s.",
+ clocks_from_ticks(getticks() - tic_total), clocks_getunit());
+ }
+
#endif /* WITH_MPI */
}
/**
- * @brief Perform a FOF search on gravity particles using the cells and applying
- * the Union-Find algorithm.
+ * @brief Compute the local size of each FOF group fragment.
*
* @param props The properties of the FOF scheme.
- * @param bh_props The properties of the black hole scheme.
- * @param constants The physical constants in internal units.
- * @param cosmo The current cosmological model.
* @param s The #space containing the particles.
- * @param dump_debug_results Are we writing txt-file debug catalogues including
- * BH-seeding info?
- * @param dump_results Do we want to write the group catalogue to a hdf5 file?
- * @param seed_black_holes Do we want to seed black holes in haloes?
*/
-void fof_search_tree(struct fof_props *props,
- const struct black_holes_props *bh_props,
- const struct phys_const *constants,
- const struct cosmology *cosmo, struct space *s,
- const int dump_results, const int dump_debug_results,
- const int seed_black_holes) {
-
- const size_t nr_gparts = s->nr_gparts;
- const size_t min_group_size = props->min_group_size;
-#ifndef WITHOUT_GROUP_PROPS
- const size_t group_id_offset = props->group_id_offset;
- const size_t group_id_default = props->group_id_default;
-#endif
+void fof_compute_local_sizes(struct fof_props *props, struct space *s) {
-#ifdef WITH_MPI
- const int nr_nodes = s->e->nr_nodes;
-#endif
- struct gpart *gparts = s->gparts;
- size_t *group_index, *group_size;
- long long num_groups = 0, num_parts_in_groups = 0, max_group_size = 0;
const int verbose = s->e->verbose;
- const ticks tic_total = getticks();
- char output_file_name[PARSER_MAX_LINE_SIZE];
- snprintf(output_file_name, PARSER_MAX_LINE_SIZE, "%s", props->base_name);
+ struct gpart *gparts = s->gparts;
+ const size_t nr_gparts = s->nr_gparts;
- if (verbose)
- message("Searching %zu gravity particles for links with l_x: %lf",
- nr_gparts, sqrt(props->l_x2));
+ const ticks tic_total = getticks();
if (engine_rank == 0 && verbose)
message("Size of hash table element: %ld", sizeof(hashmap_element_t));
#ifdef WITH_MPI
- /* Reset global variable */
- node_offset = 0;
+ const ticks comms_tic = getticks();
/* Determine number of gparts on lower numbered MPI ranks */
+ const long long nr_gparts_local = s->nr_gparts;
long long nr_gparts_cumulative;
- long long nr_gparts_local = s->nr_gparts;
-
- const ticks comms_tic = getticks();
-
MPI_Scan(&nr_gparts_local, &nr_gparts_cumulative, 1, MPI_LONG_LONG, MPI_SUM,
MPI_COMM_WORLD);
@@ -2860,13 +3767,12 @@ void fof_search_tree(struct fof_props *props,
message("MPI_Scan Imbalance took: %.3f %s.",
clocks_from_ticks(getticks() - comms_tic), clocks_getunit());
+ /* Reset global variable containing the rank particle count offset */
node_offset = nr_gparts_cumulative - nr_gparts_local;
#endif
- /* Local copy of the arrays */
- group_index = props->group_index;
- group_size = props->group_size;
-
+ /* Compute the group sizes of the local fragments
+ * (in non-MPI land that is the final group size of the haloes) */
const ticks tic_calc_group_size = getticks();
threadpool_map(&s->e->threadpool, fof_calc_group_size_mapper, gparts,
@@ -2877,31 +3783,52 @@ void fof_search_tree(struct fof_props *props,
clocks_from_ticks(getticks() - tic_calc_group_size),
clocks_getunit());
-#ifdef WITH_MPI
- if (nr_nodes > 1) {
+ if (verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic_total),
+ clocks_getunit());
+}
- const ticks tic_mpi = getticks();
+/**
+ * @brief Compute all the group properties
+ *
+ * @param props The properties of the FOF scheme.
+ * @param bh_props The properties of the black hole scheme.
+ * @param constants The physical constants in internal units.
+ * @param cosmo The current cosmological model.
+ * @param s The #space containing the particles.
+ * @param dump_debug_results Are we writing txt-file debug catalogues including
+ * BH-seeding info?
+ * @param dump_results Do we want to write the group catalogue to a hdf5 file?
+ * @param seed_black_holes Do we want to seed black holes in haloes?
+ */
+void fof_compute_group_props(struct fof_props *props,
+ const struct black_holes_props *bh_props,
+ const struct phys_const *constants,
+ const struct cosmology *cosmo, struct space *s,
+ const int dump_results,
+ const int dump_debug_results,
+ const int seed_black_holes) {
- /* Search for group links across MPI domains. */
- fof_search_foreign_cells(props, s);
+ const int verbose = s->e->verbose;
+#ifdef WITH_MPI
+ const int nr_nodes = s->e->nr_nodes;
+#endif
+ const ticks tic_total = getticks();
- if (verbose) {
- message("fof_search_foreign_cells() took (FOF SCALING): %.3f %s.",
- clocks_from_ticks(getticks() - tic_mpi), clocks_getunit());
+ struct gpart *gparts = s->gparts;
+ const size_t nr_gparts = s->nr_gparts;
- message(
- "fof_search_foreign_cells() + calc_group_size took (FOF SCALING): "
- "%.3f %s.",
- clocks_from_ticks(getticks() - tic_total), clocks_getunit());
- }
- }
-#endif
+ const size_t min_group_size = props->min_group_size;
+ const size_t group_id_offset = props->group_id_offset;
+ const size_t group_id_default = props->group_id_default;
size_t num_groups_local = 0;
-#ifndef WITHOUT_GROUP_PROPS
size_t num_parts_in_groups_local = 0;
size_t max_group_size_local = 0;
-#endif
+
+ /* Local copy of the arrays */
+ size_t *restrict group_index = props->group_index;
+ size_t *restrict group_size = props->group_size;
const ticks tic_num_groups_calc = getticks();
@@ -2917,7 +3844,6 @@ void fof_search_tree(struct fof_props *props,
num_groups_local++;
#endif
-#ifndef WITHOUT_GROUP_PROPS
/* Find the total number of particles in groups. */
if (group_size[i] >= min_group_size)
num_parts_in_groups_local += group_size[i];
@@ -2925,7 +3851,6 @@ void fof_search_tree(struct fof_props *props,
/* Find the largest group. */
if (group_size[i] > max_group_size_local)
max_group_size_local = group_size[i];
-#endif
}
if (verbose)
@@ -2934,9 +3859,8 @@ void fof_search_tree(struct fof_props *props,
"%s.",
clocks_from_ticks(getticks() - tic_num_groups_calc), clocks_getunit());
- /* Sort the groups in descending order based upon size and re-label their
- * IDs 0-num_groups. */
-#ifndef WITHOUT_GROUP_PROPS
+ /* Sort the groups in descending order based upon size and re-label their
+ * IDs 0-num_groups. */
struct group_length *high_group_sizes = NULL;
int group_count = 0;
@@ -2961,9 +3885,9 @@ void fof_search_tree(struct fof_props *props,
}
ticks tic = getticks();
-#endif /* #ifndef WITHOUT_GROUP_PROPS */
/* Find global properties. */
+ long long num_groups = 0, num_parts_in_groups = 0, max_group_size = 0;
#ifdef WITH_MPI
MPI_Allreduce(&num_groups_local, &num_groups, 1, MPI_LONG_LONG_INT, MPI_SUM,
MPI_COMM_WORLD);
@@ -2973,24 +3897,18 @@ void fof_search_tree(struct fof_props *props,
clocks_from_ticks(getticks() - tic_num_groups_calc),
clocks_getunit());
-#ifndef WITHOUT_GROUP_PROPS
MPI_Reduce(&num_parts_in_groups_local, &num_parts_in_groups, 1,
MPI_LONG_LONG_INT, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&max_group_size_local, &max_group_size, 1, MPI_LONG_LONG_INT,
MPI_MAX, 0, MPI_COMM_WORLD);
-#endif /* #ifndef WITHOUT_GROUP_PROPS */
#else
num_groups = num_groups_local;
-#ifndef WITHOUT_GROUP_PROPS
num_parts_in_groups = num_parts_in_groups_local;
max_group_size = max_group_size_local;
-#endif /* #ifndef WITHOUT_GROUP_PROPS */
#endif /* WITH_MPI */
props->num_groups = num_groups;
-#ifndef WITHOUT_GROUP_PROPS
-
/* Find number of groups on lower numbered MPI ranks */
#ifdef WITH_MPI
long long nglocal = num_groups_local;
@@ -3158,6 +4076,9 @@ void fof_search_tree(struct fof_props *props,
if (swift_memalign("fof_group_mass", (void **)&props->group_mass, 32,
num_groups_local * sizeof(double)) != 0)
error("Failed to allocate list of group masses for FOF search.");
+ if (swift_memalign("fof_group_size", (void **)&props->final_group_size, 32,
+ num_groups_local * sizeof(long long)) != 0)
+ error("Failed to allocate list of group masses for FOF search.");
if (swift_memalign("fof_group_centre_of_mass",
(void **)&props->group_centre_of_mass, 32,
num_groups_local * 3 * sizeof(double)) != 0)
@@ -3168,6 +4089,7 @@ void fof_search_tree(struct fof_props *props,
error("Failed to allocate list of group first positions for FOF search.");
bzero(props->group_mass, num_groups_local * sizeof(double));
+ bzero(props->final_group_size, num_groups_local * sizeof(long long));
bzero(props->group_centre_of_mass, num_groups_local * 3 * sizeof(double));
for (size_t i = 0; i < 3 * num_groups_local; i++) {
props->group_first_position[i] = -FLT_MAX;
@@ -3202,8 +4124,9 @@ void fof_search_tree(struct fof_props *props,
free(num_on_node);
free(first_on_node);
#else
- fof_calc_group_mass(props, s, seed_black_holes, num_groups_local, 0, NULL,
- NULL, props->group_mass);
+ fof_calc_group_mass(props, s, seed_black_holes, num_groups_local,
+ /*num_groups_prev=*/0, /*num_on_node=*/NULL,
+ /*first_on_node=*/NULL, props->group_mass);
#endif
/* Finalise the group data before dump */
@@ -3224,6 +4147,10 @@ void fof_search_tree(struct fof_props *props,
}
if (dump_debug_results) {
+
+ char output_file_name[PARSER_MAX_LINE_SIZE];
+ snprintf(output_file_name, PARSER_MAX_LINE_SIZE, "%s", props->base_name);
+
#ifdef WITH_MPI
snprintf(output_file_name + strlen(output_file_name), FILENAME_BUFFER_SIZE,
"_mpi.dat");
@@ -3244,17 +4171,21 @@ void fof_search_tree(struct fof_props *props,
/* Free the left-overs */
swift_free("fof_high_group_sizes", high_group_sizes);
swift_free("fof_group_mass", props->group_mass);
+ swift_free("fof_group_size", props->final_group_size);
swift_free("fof_group_centre_of_mass", props->group_centre_of_mass);
swift_free("fof_group_first_position", props->group_first_position);
swift_free("fof_max_part_density_index", props->max_part_density_index);
swift_free("fof_max_part_density", props->max_part_density);
props->group_mass = NULL;
+ props->final_group_size = NULL;
props->group_centre_of_mass = NULL;
props->max_part_density_index = NULL;
props->max_part_density = NULL;
-#endif /* #ifndef WITHOUT_GROUP_PROPS */
+ swift_free("fof_distance", props->distance_to_link);
swift_free("fof_group_index", props->group_index);
+ swift_free("fof_attach_index", props->attach_index);
+ swift_free("fof_found_attach", props->found_attachable_link);
swift_free("fof_group_size", props->group_size);
props->group_index = NULL;
props->group_size = NULL;
@@ -3286,6 +4217,7 @@ void fof_struct_dump(const struct fof_props *props, FILE *stream) {
temp.group_index = NULL;
temp.group_size = NULL;
temp.group_mass = NULL;
+ temp.final_group_size = NULL;
temp.group_centre_of_mass = NULL;
temp.max_part_density_index = NULL;
temp.max_part_density = NULL;
diff --git a/src/fof.h b/src/fof.h
index bff0c0a3d0f0933866b7d3949ebf5b5f6bdc8b30..c7213e8532c7ff614c45d18ef8ea1ee8df1a5e7c 100644
--- a/src/fof.h
+++ b/src/fof.h
@@ -25,6 +25,7 @@
/* Local headers */
#include "align.h"
#include "parser.h"
+#include "part_type.h"
/* Avoid cyclic inclusions */
struct cell;
@@ -68,6 +69,12 @@ struct fof_props {
/*! The base name of the output file */
char base_name[PARSER_MAX_LINE_SIZE];
+ /*! The types of particles to use for linking */
+ int fof_linking_types[swift_type_count];
+
+ /*! The types of particles to use for attaching */
+ int fof_attach_types[swift_type_count];
+
/* ------------ Group properties ----------------- */
/*! Number of groups */
@@ -80,9 +87,21 @@ struct fof_props {
/*! Index of the root particle of the group a given gpart belongs to. */
size_t *group_index;
+ /*! Index of the root particle of the group a given gpart is attached to. */
+ size_t *attach_index;
+
+ /*! Has the particle found a linkable to attach to? */
+ char *found_attachable_link;
+
+ /*! For attachable particles: distance to the current nearest linkable part */
+ float *distance_to_link;
+
/*! Size of the group a given gpart belongs to. */
size_t *group_size;
+ /*! Final size of the group a given gpart belongs to. */
+ long long *final_group_size;
+
/*! Mass of the group a given gpart belongs to. */
double *group_mass;
@@ -147,6 +166,7 @@ struct fof_final_index {
struct fof_final_mass {
size_t global_root;
double group_mass;
+ long long final_group_size;
double first_position[3];
double centre_of_mass[3];
long long max_part_density_index;
@@ -171,14 +191,20 @@ void fof_init(struct fof_props *props, struct swift_params *params,
const struct phys_const *phys_const, const struct unit_system *us,
const int stand_alone_fof);
void fof_create_mpi_types(void);
-void fof_allocate(const struct space *s, const long long total_nr_DM_particles,
- struct fof_props *props);
-void fof_search_tree(struct fof_props *props,
- const struct black_holes_props *bh_props,
- const struct phys_const *constants,
- const struct cosmology *cosmo, struct space *s,
- const int dump_results, const int dump_debug_results,
- const int seed_black_holes);
+void fof_allocate(const struct space *s, struct fof_props *props);
+void fof_compute_local_sizes(struct fof_props *props, struct space *s);
+void fof_search_foreign_cells(struct fof_props *props, const struct space *s);
+void fof_link_attachable_particles(struct fof_props *props,
+ const struct space *s);
+void fof_finalise_attachables(struct fof_props *props, const struct space *s);
+void fof_link_foreign_fragments(struct fof_props *props, const struct space *s);
+void fof_compute_group_props(struct fof_props *props,
+ const struct black_holes_props *bh_props,
+ const struct phys_const *constants,
+ const struct cosmology *cosmo, struct space *s,
+ const int dump_results,
+ const int dump_debug_results,
+ const int seed_black_holes);
void rec_fof_search_self(const struct fof_props *props, const double dim[3],
const double search_r2, const int periodic,
const struct gpart *const space_gparts,
@@ -187,6 +213,16 @@ void rec_fof_search_pair(const struct fof_props *props, const double dim[3],
const double search_r2, const int periodic,
const struct gpart *const space_gparts,
struct cell *restrict ci, struct cell *restrict cj);
+void rec_fof_attach_self(const struct fof_props *props, const double dim[3],
+ const double search_r2, const int periodic,
+ const struct gpart *const space_gparts,
+ const size_t nr_gparts, struct cell *c);
+void rec_fof_attach_pair(const struct fof_props *props, const double dim[3],
+ const double search_r2, const int periodic,
+ const struct gpart *const space_gparts,
+ const size_t nr_gparts, struct cell *restrict ci,
+ struct cell *restrict cj, const int ci_local,
+ const int cj_local);
void fof_struct_dump(const struct fof_props *props, FILE *stream);
void fof_struct_restore(struct fof_props *props, FILE *stream);
diff --git a/src/fof_catalogue_io.c b/src/fof_catalogue_io.c
index 12308fb75102469890ba095b30627e1ab85c9026..8c1b20d4bf4c96cfb26def11608c9627bfd2495b 100644
--- a/src/fof_catalogue_io.c
+++ b/src/fof_catalogue_io.c
@@ -35,7 +35,7 @@
void write_fof_hdf5_header(hid_t h_file, const struct engine* e,
const long long num_groups_total,
const long long num_groups_this_file,
- const struct fof_props* props) {
+ const struct fof_props* props, const int virtual) {
/* Open header to write simulation properties */
hid_t h_grp =
@@ -117,7 +117,7 @@ void write_fof_hdf5_header(hid_t h_file, const struct engine* e,
io_write_attribute_i(h_grp, "NumFilesPerSnapshot", e->nr_nodes);
io_write_attribute_i(h_grp, "ThisFile", e->nodeID);
io_write_attribute_s(h_grp, "SelectOutput", "Default");
- io_write_attribute_i(h_grp, "Virtual", 0);
+ io_write_attribute_i(h_grp, "Virtual", virtual);
const int to_write[swift_type_count] = {0};
io_write_attribute(h_grp, "CanHaveTypes", INT, to_write, swift_type_count);
io_write_attribute_s(h_grp, "OutputType", "FOF");
@@ -133,7 +133,221 @@ void write_fof_hdf5_header(hid_t h_file, const struct engine* e,
ic_info_write_hdf5(e->ics_metadata, h_file);
/* Write all the meta-data */
- io_write_meta_data(h_file, e, e->internal_units, e->snapshot_units);
+ io_write_meta_data(h_file, e, e->internal_units, e->snapshot_units,
+ /*fof=*/1);
+}
+
+void write_virtual_fof_hdf5_array(
+ const struct engine* e, hid_t grp, const char* fileName_base,
+ const char* partTypeGroupName, const struct io_props props,
+ const size_t N_total, const long long* N_counts,
+ const enum lossy_compression_schemes lossy_compression,
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units) {
+
+#if H5_VERSION_GE(1, 10, 0)
+
+ /* Create data space */
+ hid_t h_space;
+ if (N_total > 0)
+ h_space = H5Screate(H5S_SIMPLE);
+ else
+ h_space = H5Screate(H5S_NULL);
+
+ if (h_space < 0)
+ error("Error while creating data space for field '%s'.", props.name);
+
+ int rank = 0;
+ hsize_t shape[2];
+ hsize_t source_shape[2];
+ hsize_t start[2] = {0, 0};
+ hsize_t count[2];
+ if (props.dimension > 1) {
+ rank = 2;
+ shape[0] = N_total;
+ shape[1] = props.dimension;
+ source_shape[0] = 0;
+ source_shape[1] = props.dimension;
+ count[0] = 0;
+ count[1] = props.dimension;
+
+ } else {
+ rank = 1;
+ shape[0] = N_total;
+ shape[1] = 0;
+ source_shape[0] = 0;
+ source_shape[1] = 0;
+ count[0] = 0;
+ count[1] = 0;
+ }
+
+ /* Change shape of data space */
+ hid_t h_err = H5Sset_extent_simple(h_space, rank, shape, NULL);
+ if (h_err < 0)
+ error("Error while changing data space shape for field '%s'.", props.name);
+
+ /* Dataset type */
+ hid_t h_type = H5Tcopy(io_hdf5_type(props.type));
+
+ /* Dataset properties */
+ hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);
+
+ /* Create filters and set compression level if we have something to write */
+ char comp_buffer[32] = "None";
+
+ /* The name of the dataset to map to in the other files */
+ char source_dataset_name[256];
+ sprintf(source_dataset_name, "Groups/%s", props.name);
+
+ /* Construct a relative base name */
+ char fileName_relative_base[256];
+ int pos_last_slash = strlen(fileName_base) - 1;
+ for (/* */; pos_last_slash >= 0; --pos_last_slash)
+ if (fileName_base[pos_last_slash] == '/') break;
+
+ sprintf(fileName_relative_base, "%s", &fileName_base[pos_last_slash + 1]);
+
+ /* Create all the virtual mappings */
+ for (int i = 0; i < e->nr_nodes; ++i) {
+
+ /* Get the number of particles of this type written on this rank */
+ count[0] = N_counts[i];
+
+ /* Select the space in the virtual file */
+ h_err = H5Sselect_hyperslab(h_space, H5S_SELECT_SET, start, /*stride=*/NULL,
+ count, /*block=*/NULL);
+ if (h_err < 0) error("Error selecting hyper-slab in the virtual file");
+
+ /* Select the space in the (already existing) source file */
+ source_shape[0] = count[0];
+ hid_t h_source_space = H5Screate_simple(rank, source_shape, NULL);
+ if (h_source_space < 0) error("Error creating space in the source file");
+
+ char fileName[1024];
+ sprintf(fileName, "%s.%d.hdf5", fileName_relative_base, i);
+
+ /* Make the virtual link */
+ h_err = H5Pset_virtual(h_prop, h_space, fileName, source_dataset_name,
+ h_source_space);
+ if (h_err < 0) error("Error setting the virtual properties");
+
+ H5Sclose(h_source_space);
+
+ /* Move to the next slab (i.e. next file) */
+ start[0] += count[0];
+ }
+
+ /* Create virtual dataset */
+ const hid_t h_data = H5Dcreate(grp, props.name, h_type, h_space, H5P_DEFAULT,
+ h_prop, H5P_DEFAULT);
+ if (h_data < 0) error("Error while creating dataspace '%s'.", props.name);
+
+ /* Write unit conversion factors for this data set */
+ char buffer[FIELD_BUFFER_SIZE] = {0};
+ units_cgs_conversion_string(buffer, snapshot_units, props.units,
+ props.scale_factor_exponent);
+ float baseUnitsExp[5];
+ units_get_base_unit_exponents_array(baseUnitsExp, props.units);
+ io_write_attribute_f(h_data, "U_M exponent", baseUnitsExp[UNIT_MASS]);
+ io_write_attribute_f(h_data, "U_L exponent", baseUnitsExp[UNIT_LENGTH]);
+ io_write_attribute_f(h_data, "U_t exponent", baseUnitsExp[UNIT_TIME]);
+ io_write_attribute_f(h_data, "U_I exponent", baseUnitsExp[UNIT_CURRENT]);
+ io_write_attribute_f(h_data, "U_T exponent", baseUnitsExp[UNIT_TEMPERATURE]);
+ io_write_attribute_f(h_data, "h-scale exponent", 0.f);
+ io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
+ io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
+ io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
+
+ /* Write the actual number this conversion factor corresponds to */
+ const double factor =
+ units_cgs_conversion_factor(snapshot_units, props.units);
+ io_write_attribute_d(
+ h_data,
+ "Conversion factor to CGS (not including cosmological corrections)",
+ factor);
+ io_write_attribute_d(
+ h_data,
+ "Conversion factor to physical CGS (including cosmological corrections)",
+ factor * pow(e->cosmology->a, props.scale_factor_exponent));
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (strlen(props.description) == 0)
+ error("Invalid (empty) description of the field '%s'", props.name);
+#endif
+
+ /* Write the full description */
+ io_write_attribute_s(h_data, "Description", props.description);
+
+ /* Close everything */
+ H5Tclose(h_type);
+ H5Pclose(h_prop);
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+
+#endif
+}
+
+void write_fof_virtual_file(const struct fof_props* props,
+ const size_t num_groups_total,
+ const long long* N_counts, const struct engine* e) {
+#if H5_VERSION_GE(1, 10, 0)
+
+ /* Create the filename */
+ char file_name[512];
+ char file_name_base[512];
+ char subdir_name[265];
+ sprintf(subdir_name, "%s_%04i", props->base_name, e->snapshot_output_count);
+ const char* base = basename(subdir_name);
+ sprintf(file_name_base, "%s/%s", subdir_name, base);
+ sprintf(file_name, "%s/%s.hdf5", subdir_name, base);
+
+ /* Open HDF5 file with the chosen parameters */
+ hid_t h_file = H5Fcreate(file_name, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_file < 0) error("Error while opening file '%s'.", file_name);
+
+ /* Start by writing the header */
+ write_fof_hdf5_header(h_file, e, num_groups_total, num_groups_total, props,
+ /*virtual=*/1);
+ hid_t h_grp =
+ H5Gcreate(h_file, "/Groups", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating groups group.\n");
+
+ struct io_props output_prop;
+ output_prop = io_make_output_field_("Masses", DOUBLE, 1, UNIT_CONV_MASS, 0.f,
+ (char*)props->group_mass, sizeof(double),
+ "FOF group masses");
+ write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
+ num_groups_total, N_counts,
+ compression_write_lossless, e->internal_units,
+ e->snapshot_units);
+ output_prop =
+ io_make_output_field_("Centres", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f,
+ (char*)props->group_centre_of_mass,
+ 3 * sizeof(double), "FOF group centres of mass");
+ write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
+ num_groups_total, N_counts,
+ compression_write_lossless, e->internal_units,
+ e->snapshot_units);
+ output_prop = io_make_output_field_(
+ "GroupIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
+ (char*)props->group_index, sizeof(size_t), "FOF group IDs");
+ write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
+ num_groups_total, N_counts,
+ compression_write_lossless, e->internal_units,
+ e->snapshot_units);
+ output_prop =
+ io_make_output_field_("Sizes", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
+ (char*)props->final_group_size, sizeof(long long),
+ "FOF group length (number of particles)");
+ write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
+ num_groups_total, N_counts,
+ compression_write_lossless, e->internal_units,
+ e->snapshot_units);
+
+ /* Close everything */
+ H5Gclose(h_grp);
+ H5Fclose(h_file);
+#endif
}
void write_fof_hdf5_array(
@@ -316,10 +530,16 @@ void write_fof_hdf5_catalogue(const struct fof_props* props,
#ifdef WITH_MPI
MPI_Allreduce(&num_groups, &num_groups_total, 1, MPI_LONG_LONG, MPI_SUM,
MPI_COMM_WORLD);
+
+ /* Rank 0 collects the number of groups written by each rank */
+ long long* N_counts = (long long*)malloc(e->nr_nodes * sizeof(long long));
+ MPI_Gather(&num_groups_local, 1, MPI_LONG_LONG_INT, N_counts, 1,
+ MPI_LONG_LONG_INT, 0, MPI_COMM_WORLD);
#endif
/* Start by writing the header */
- write_fof_hdf5_header(h_file, e, num_groups_total, num_groups_local, props);
+ write_fof_hdf5_header(h_file, e, num_groups_total, num_groups_local, props,
+ /*virtual=*/0);
hid_t h_grp =
H5Gcreate(h_file, "/Groups", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
@@ -345,9 +565,10 @@ void write_fof_hdf5_catalogue(const struct fof_props* props,
write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop,
num_groups_local, compression_write_lossless,
e->internal_units, e->snapshot_units);
- output_prop = io_make_output_field_(
- "Sizes", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, (char*)props->group_size,
- sizeof(size_t), "FOF group length (number of particles)");
+ output_prop =
+ io_make_output_field_("Sizes", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
+ (char*)props->final_group_size, sizeof(long long),
+ "FOF group length (number of particles)");
write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop,
num_groups_local, compression_write_lossless,
e->internal_units, e->snapshot_units);
@@ -357,6 +578,16 @@ void write_fof_hdf5_catalogue(const struct fof_props* props,
H5Fclose(h_file);
#ifdef WITH_MPI
+#if H5_VERSION_GE(1, 10, 0)
+
+ /* Write the virtual meta-file */
+ if (e->nodeID == 0)
+ write_fof_virtual_file(props, num_groups_total, N_counts, e);
+#endif
+
+ /* Free the counts-per-rank array */
+ free(N_counts);
+
MPI_Barrier(MPI_COMM_WORLD);
#endif
}
diff --git a/src/hashmap.h b/src/hashmap.h
index 853f3965315b0a079799389a93e3869725ee35b7..6381e26e036463fbb672dd0552656bf6fe7fa44a 100644
--- a/src/hashmap.h
+++ b/src/hashmap.h
@@ -50,6 +50,7 @@ typedef size_t hashmap_mask_t;
#ifndef hashmap_value_t
typedef struct _hashmap_struct {
long long value_st;
+ long long value_ll;
float value_flt;
double value_dbl;
double value_array_dbl[3];
diff --git a/src/line_of_sight.c b/src/line_of_sight.c
index 01680ab80721e183b33c9fad17ec1968fe494eec..40e71f960224c60329f68e04f206154f71c29d73 100644
--- a/src/line_of_sight.c
+++ b/src/line_of_sight.c
@@ -551,7 +551,8 @@ void write_hdf5_header(hid_t h_file, const struct engine *e,
ic_info_write_hdf5(e->ics_metadata, h_file);
/* Write all the meta-data */
- io_write_meta_data(h_file, e, e->internal_units, e->snapshot_units);
+ io_write_meta_data(h_file, e, e->internal_units, e->snapshot_units,
+ /*fof=*/0);
/* Print the LOS properties */
h_grp = H5Gcreate(h_file, "/LineOfSightParameters", H5P_DEFAULT, H5P_DEFAULT,
diff --git a/src/parallel_io.c b/src/parallel_io.c
index bc81e1bbc433eb749d28858a290bbe12a2748b0f..97564cecd99038d513ee8bc0e9279f4dfa974ea8 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1148,6 +1148,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
* @param numFields The number of fields to write for each particle type.
* @param internal_units The #unit_system used internally.
* @param snapshot_units The #unit_system used in the snapshots.
+ * @param fof Is this a snapshot related to a stand-alone FOF call?
* @param subsample_any Are any fields being subsampled?
* @param subsample_fraction The subsampling fraction of each particle type.
*/
@@ -1158,7 +1159,7 @@ void prepare_file(struct engine* e, const char* fileName,
const int numFields[swift_type_count],
const char current_selection_name[FIELD_BUFFER_SIZE],
const struct unit_system* internal_units,
- const struct unit_system* snapshot_units,
+ const struct unit_system* snapshot_units, const int fof,
const int subsample_any,
const float subsample_fraction[swift_type_count]) {
@@ -1294,7 +1295,7 @@ void prepare_file(struct engine* e, const char* fileName,
ic_info_write_hdf5(e->ics_metadata, h_file);
/* Write all the meta-data */
- io_write_meta_data(h_file, e, internal_units, snapshot_units);
+ io_write_meta_data(h_file, e, internal_units, snapshot_units, fof);
/* Loop over all particle types */
for (int ptype = 0; ptype < swift_type_count; ptype++) {
@@ -1423,6 +1424,7 @@ void prepare_file(struct engine* e, const char* fileName,
* @param e The engine containing all the system.
* @param internal_units The #unit_system used internally
* @param snapshot_units The #unit_system used in the snapshots
+ * @param fof Is this a snapshot related to a stand-alone FOF call?
* @param mpi_rank The MPI rank of this node.
* @param mpi_size The number of MPI ranks.
* @param comm The MPI communicator.
@@ -1439,8 +1441,8 @@ void prepare_file(struct engine* e, const char* fileName,
void write_output_parallel(struct engine* e,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units,
- const int mpi_rank, const int mpi_size,
- MPI_Comm comm, MPI_Info info) {
+ const int fof, const int mpi_rank,
+ const int mpi_size, MPI_Comm comm, MPI_Info info) {
const struct part* parts = e->s->parts;
const struct xpart* xparts = e->s->xparts;
@@ -1624,7 +1626,7 @@ void write_output_parallel(struct engine* e,
/* Rank 0 prepares the file */
if (mpi_rank == 0)
prepare_file(e, fileName, xmfFileName, N_total, to_write, numFields,
- current_selection_name, internal_units, snapshot_units,
+ current_selection_name, internal_units, snapshot_units, fof,
subsample_any, subsample_fraction);
MPI_Barrier(MPI_COMM_WORLD);
diff --git a/src/parallel_io.h b/src/parallel_io.h
index aaac00aa4c8b1153fe844e61b9146843834030cf..98be1278f5889a2b3f61459ba0ef5fe248a12860 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -52,8 +52,8 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
void write_output_parallel(struct engine* e,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units,
- int mpi_rank, int mpi_size, MPI_Comm comm,
- MPI_Info info);
+ const int fof, int mpi_rank, int mpi_size,
+ MPI_Comm comm, MPI_Info info);
#endif
#endif /* SWIFT_PARALLEL_IO_H */
diff --git a/src/runner.h b/src/runner.h
index 4ccd3cec4ac7e6ef1edffb0e81ec16da54d15145..8452b41879b0f6073091267a2205197d1978fbfe 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -124,9 +124,12 @@ void runner_do_grav_mesh(struct runner *r, struct cell *c, int timer);
void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
void runner_do_grav_fft(struct runner *r, int timer);
void runner_do_csds(struct runner *r, struct cell *c, int timer);
-void runner_do_fof_self(struct runner *r, struct cell *c, int timer);
-void runner_do_fof_pair(struct runner *r, struct cell *ci, struct cell *cj,
- int timer);
+void runner_do_fof_search_self(struct runner *r, struct cell *c, int timer);
+void runner_do_fof_search_pair(struct runner *r, struct cell *ci,
+ struct cell *cj, int timer);
+void runner_do_fof_attach_self(struct runner *r, struct cell *c, int timer);
+void runner_do_fof_attach_pair(struct runner *r, struct cell *ci,
+ struct cell *cj, int timer);
void runner_do_rt_ghost1(struct runner *r, struct cell *c, int timer);
void runner_do_rt_ghost2(struct runner *r, struct cell *c, int timer);
void runner_do_rt_tchem(struct runner *r, struct cell *c, int timer);
diff --git a/src/runner_main.c b/src/runner_main.c
index 8c32c7d709430ecbbfbe64856db4e287b0620849..66a0074b9d917c6c8c2f455ad68aef18c1853d31 100644
--- a/src/runner_main.c
+++ b/src/runner_main.c
@@ -568,10 +568,16 @@ void *runner_main(void *data) {
runner_do_sink_formation(r, t->ci);
break;
case task_type_fof_self:
- runner_do_fof_self(r, t->ci, 1);
+ runner_do_fof_search_self(r, t->ci, 1);
break;
case task_type_fof_pair:
- runner_do_fof_pair(r, t->ci, t->cj, 1);
+ runner_do_fof_search_pair(r, t->ci, t->cj, 1);
+ break;
+ case task_type_fof_attach_self:
+ runner_do_fof_attach_self(r, t->ci, 1);
+ break;
+ case task_type_fof_attach_pair:
+ runner_do_fof_attach_pair(r, t->ci, t->cj, 1);
break;
case task_type_neutrino_weight:
runner_do_neutrino_weighting(r, ci, 1);
diff --git a/src/runner_others.c b/src/runner_others.c
index c70bc0ccb511b5ad39eb2124e764c2332aeb476c..a79f2ce06563e830e4e4a8d8018763cfbb541dbd 100644
--- a/src/runner_others.c
+++ b/src/runner_others.c
@@ -956,7 +956,7 @@ void runner_do_csds(struct runner *r, struct cell *c, int timer) {
* @param c cell
* @param timer 1 if the time is to be recorded.
*/
-void runner_do_fof_self(struct runner *r, struct cell *c, int timer) {
+void runner_do_fof_search_self(struct runner *r, struct cell *c, int timer) {
#ifdef WITH_FOF
@@ -986,8 +986,8 @@ void runner_do_fof_self(struct runner *r, struct cell *c, int timer) {
* @param cj cell j
* @param timer 1 if the time is to be recorded.
*/
-void runner_do_fof_pair(struct runner *r, struct cell *ci, struct cell *cj,
- int timer) {
+void runner_do_fof_search_pair(struct runner *r, struct cell *ci,
+ struct cell *cj, int timer) {
#ifdef WITH_FOF
@@ -1013,6 +1013,68 @@ void runner_do_fof_pair(struct runner *r, struct cell *ci, struct cell *cj,
#endif
}
+/**
+ * @brief Recursively search for FOF groups in a single cell.
+ *
+ * @param r runner task
+ * @param c cell
+ * @param timer 1 if the time is to be recorded.
+ */
+void runner_do_fof_attach_self(struct runner *r, struct cell *c, int timer) {
+
+#ifdef WITH_FOF
+
+ TIMER_TIC;
+
+ const struct engine *e = r->e;
+ struct space *s = e->s;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
+ const int periodic = s->periodic;
+ const struct gpart *const gparts = s->gparts;
+ const double attach_r2 = e->fof_properties->l_x2;
+
+ rec_fof_attach_self(e->fof_properties, dim, attach_r2, periodic, gparts,
+ s->nr_gparts, c);
+
+ if (timer) TIMER_TOC(timer_fof_self);
+
+#else
+ error("SWIFT was not compiled with FOF enabled!");
+#endif
+}
+
+/**
+ * @brief Recursively search for FOF groups between a pair of cells.
+ *
+ * @param r runner task
+ * @param ci cell i
+ * @param cj cell j
+ * @param timer 1 if the time is to be recorded.
+ */
+void runner_do_fof_attach_pair(struct runner *r, struct cell *ci,
+ struct cell *cj, int timer) {
+
+#ifdef WITH_FOF
+
+ TIMER_TIC;
+
+ const struct engine *e = r->e;
+ struct space *s = e->s;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
+ const int periodic = s->periodic;
+ const struct gpart *const gparts = s->gparts;
+ const double attach_r2 = e->fof_properties->l_x2;
+
+ rec_fof_attach_pair(e->fof_properties, dim, attach_r2, periodic, gparts,
+ s->nr_gparts, ci, cj, e->nodeID == ci->nodeID,
+ e->nodeID == cj->nodeID);
+
+ if (timer) TIMER_TOC(timer_fof_pair);
+#else
+ error("SWIFT was not compiled with FOF enabled!");
+#endif
+}
+
/**
* @brief Finish up the transport step and do the thermochemistry
* for radiative transfer
diff --git a/src/serial_io.c b/src/serial_io.c
index 2e39dc4d2dc53707cb2a792c945d3b579a04fc51..586c188fea37f7d394e592cdb2ef7e0e485cda1c 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -945,6 +945,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
* @param e The engine containing all the system.
* @param internal_units The #unit_system used internally
* @param snapshot_units The #unit_system used in the snapshots
+ * @param fof Is this a snapshot related to a stand-alone FOF call?
* @param mpi_rank The MPI rank of this node.
* @param mpi_size The number of MPI ranks.
* @param comm The MPI communicator.
@@ -961,8 +962,8 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
void write_output_serial(struct engine* e,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units,
- const int mpi_rank, const int mpi_size, MPI_Comm comm,
- MPI_Info info) {
+ const int fof, const int mpi_rank, const int mpi_size,
+ MPI_Comm comm, MPI_Info info) {
hid_t h_file = 0, h_grp = 0;
int numFiles = 1;
@@ -1263,7 +1264,7 @@ void write_output_serial(struct engine* e,
ic_info_write_hdf5(e->ics_metadata, h_file);
/* Write all the meta-data */
- io_write_meta_data(h_file, e, internal_units, snapshot_units);
+ io_write_meta_data(h_file, e, internal_units, snapshot_units, fof);
/* Loop over all particle types */
for (int ptype = 0; ptype < swift_type_count; ptype++) {
diff --git a/src/serial_io.h b/src/serial_io.h
index 350da844aa034575c04f7efa432e84abb933c5ef..2a323f543f53bad03e3e661d6d8189b69d907f68 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -54,8 +54,8 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
void write_output_serial(struct engine* e,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units,
- const int mpi_rank, const int mpi_size, MPI_Comm comm,
- MPI_Info info);
+ const int fof, const int mpi_rank, const int mpi_size,
+ MPI_Comm comm, MPI_Info info);
#endif /* HAVE_HDF5 && WITH_MPI && !HAVE_PARALLEL_HDF5 */
diff --git a/src/single_io.c b/src/single_io.c
index 4fffaee3ab77b9cb5d2f65333e3b41ae7ae95b33..e52933e7f7b876a9e3b3c919871c82e36d2d87c5 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -799,6 +799,7 @@ void read_ic_single(
* @param e The engine containing all the system.
* @param internal_units The #unit_system used internally
* @param snapshot_units The #unit_system used in the snapshots
+ * @param fof Is this a snapshot related to a stand-alone FOF call?
*
* Creates an HDF5 output file and writes the particles contained
* in the engine. If such a file already exists, it is erased and replaced
@@ -810,7 +811,8 @@ void read_ic_single(
*/
void write_output_single(struct engine* e,
const struct unit_system* internal_units,
- const struct unit_system* snapshot_units) {
+ const struct unit_system* snapshot_units,
+ const int fof) {
hid_t h_file = 0, h_grp = 0;
int numFiles = 1;
@@ -1093,7 +1095,7 @@ void write_output_single(struct engine* e,
ic_info_write_hdf5(e->ics_metadata, h_file);
/* Write all the meta-data */
- io_write_meta_data(h_file, e, internal_units, snapshot_units);
+ io_write_meta_data(h_file, e, internal_units, snapshot_units, fof);
/* Now write the top-level cell structure */
long long global_offsets[swift_type_count] = {0};
diff --git a/src/single_io.h b/src/single_io.h
index 13a1e1be9c292aafc3fd8c0c5754fa32b7e53512..7d473d863b3ed0b8e73b8af1f437840adf982741 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -45,7 +45,8 @@ void read_ic_single(
void write_output_single(struct engine* e,
const struct unit_system* internal_units,
- const struct unit_system* snapshot_units);
+ const struct unit_system* snapshot_units,
+ const int fof);
#endif /* HAVE_HDF5 && !WITH_MPI */
diff --git a/src/task.c b/src/task.c
index 0fff5fd03b88f9ed86ca8ae00ee0e474ad56627e..3b6ae49f48c60900e2e983a077c393fcda844cef 100644
--- a/src/task.c
+++ b/src/task.c
@@ -110,6 +110,8 @@ const char *taskID_names[task_type_count] = {
"bh_swallow_ghost3",
"fof_self",
"fof_pair",
+ "fof_attach_self",
+ "fof_attach_pair",
"neutrino_weight",
"sink_in",
"sink_ghost1",
@@ -320,6 +322,8 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_type_csds:
case task_type_fof_self:
case task_type_fof_pair:
+ case task_type_fof_attach_self:
+ case task_type_fof_attach_pair:
case task_type_timestep:
case task_type_timestep_limiter:
case task_type_timestep_sync:
@@ -679,6 +683,17 @@ void task_unlock(struct task *t) {
#endif
break;
+ case task_type_fof_self:
+ case task_type_fof_attach_self:
+ cell_gunlocktree(ci);
+ break;
+
+ case task_type_fof_pair:
+ case task_type_fof_attach_pair:
+ cell_gunlocktree(ci);
+ cell_gunlocktree(cj);
+ break;
+
case task_type_star_formation:
cell_unlocktree(ci);
cell_sunlocktree(ci);
@@ -1072,6 +1087,24 @@ int task_lock(struct task *t) {
#endif
break;
+ case task_type_fof_self:
+ case task_type_fof_attach_self:
+ /* Lock the gpart as this this what we act on */
+ if (ci->grav.phold) return 0;
+ if (cell_glocktree(ci) != 0) return 0;
+ break;
+
+ case task_type_fof_pair:
+ case task_type_fof_attach_pair:
+ /* Lock the gpart as this this what we act on */
+ if (ci->grav.phold || cj->grav.phold) return 0;
+ if (cell_glocktree(ci) != 0) return 0;
+ if (cell_glocktree(cj) != 0) {
+ cell_gunlocktree(ci);
+ return 0;
+ }
+ break;
+
case task_type_star_formation:
/* Lock the gas, gravity and star particles */
if (ci->hydro.hold || ci->stars.hold || ci->grav.phold) return 0;
@@ -1776,6 +1809,8 @@ enum task_categories task_get_category(const struct task *t) {
case task_type_fof_self:
case task_type_fof_pair:
+ case task_type_fof_attach_self:
+ case task_type_fof_attach_pair:
return task_category_fof;
case task_type_rt_in:
diff --git a/src/task.h b/src/task.h
index 54f968b7d9c5e8f1e95ed8ef35332e44a7bb5ab9..47649d621b23f0788bf35cd540faaadc287bfe86 100644
--- a/src/task.h
+++ b/src/task.h
@@ -103,6 +103,8 @@ enum task_types {
task_type_bh_swallow_ghost3, /* Implicit */
task_type_fof_self,
task_type_fof_pair,
+ task_type_fof_attach_self,
+ task_type_fof_attach_pair,
task_type_neutrino_weight,
task_type_sink_in, /* Implicit */
task_type_sink_ghost1, /* Implicit */
diff --git a/swift.c b/swift.c
index 8c54e51d1ec2741c660e87abe96bf7211b9ac32d..1851ad38684f17ebe7b80d7c5124c7b0a6f8b6de 100644
--- a/swift.c
+++ b/swift.c
@@ -1633,7 +1633,7 @@ int main(int argc, char *argv[]) {
if (with_power)
calc_all_power_spectra(e.power_data, e.s, &e.threadpool, e.verbose);
- engine_dump_snapshot(&e);
+ engine_dump_snapshot(&e, /*fof=*/0);
}
/* Dump initial state statistics, if not working with an output list */
@@ -1876,7 +1876,7 @@ int main(int argc, char *argv[]) {
!e.stf_this_timestep)
velociraptor_invoke(&e, /*linked_with_snap=*/1);
#endif
- engine_dump_snapshot(&e);
+ engine_dump_snapshot(&e, /*fof=*/0);
#ifdef HAVE_VELOCIRAPTOR
if (with_structure_finding && e.snapshot_invoke_stf &&
e.s->gpart_group_data)
diff --git a/swift_fof.c b/swift_fof.c
index af7aefd7848ee5d4c8bfa50f67708f2e900efb1e..0930426f48a3f49a40719a55158b78d470ec0c91 100644
--- a/swift_fof.c
+++ b/swift_fof.c
@@ -716,7 +716,7 @@ int main(int argc, char *argv[]) {
if (with_sinks) e.policy |= engine_policy_sinks;
/* Write output. */
- engine_dump_snapshot(&e);
+ engine_dump_snapshot(&e, /*fof=*/1);
#ifdef WITH_MPI
MPI_Barrier(MPI_COMM_WORLD);
@@ -769,11 +769,6 @@ int main(int argc, char *argv[]) {
}
#endif
-#ifdef WITH_MPI
- if ((res = MPI_Finalize()) != MPI_SUCCESS)
- error("call to MPI_Finalize failed with error %i.", res);
-#endif
-
/* Clean everything */
cosmology_clean(&cosmo);
pm_mesh_clean(&mesh);
@@ -781,6 +776,11 @@ int main(int argc, char *argv[]) {
free(params);
free(output_options);
+#ifdef WITH_MPI
+ if ((res = MPI_Finalize()) != MPI_SUCCESS)
+ error("call to MPI_Finalize failed with error %i.", res);
+#endif
+
/* Say goodbye. */
if (myrank == 0) message("done. Bye.");
diff --git a/tests/testSelectOutput.c b/tests/testSelectOutput.c
index b6a73c54828be921a2907ab9ebd21488f4f245da..527af39edfdb43cce064c36499e6d5bb6f210c28 100644
--- a/tests/testSelectOutput.c
+++ b/tests/testSelectOutput.c
@@ -172,7 +172,7 @@ int main(int argc, char *argv[]) {
/* write output file */
message("Writing output.");
- write_output_single(&e, &us, &us);
+ write_output_single(&e, &us, &us, /*fof=*/0);
/* Clean-up */
message("Cleaning memory.");
diff --git a/tools/task_plots/swift_hardcoded_data.py b/tools/task_plots/swift_hardcoded_data.py
index 76277d49b2caffc76ebf26a5aa1f8dd2cbb382e1..cafc2233b2259be246c4acd8ad08991d94dfe34a 100644
--- a/tools/task_plots/swift_hardcoded_data.py
+++ b/tools/task_plots/swift_hardcoded_data.py
@@ -70,6 +70,8 @@ TASKTYPES = [
"bh_swallow_ghost3",
"fof_self",
"fof_pair",
+ "fof_attach_self",
+ "fof_attach_pair",
"neutrino_weight",
"sink_in",
"sink_ghost1",