swiftsim/soap

Spherical Overdensity Aperture Processor: MPI parallel Python code to compute properties of halos in SWIFT n-body simulations

SWIFT: SPH WIth Fine-grained inter-dependent Tasking

SWIFT is a gravity and SPH solver designed to run cosmological simulations on peta-scale machines, scaling well up to 10's of thousands of compute node.

More general information about SWIFT is available on the project webpages.

For information on how to run SWIFT, please consult the onboarding guide available here. This includes dependencies, and a few examples to get you going.

We suggest that you use the latest release branch of SWIFT, rather than the current master branch as this will change rapidly. We do, however, like to ensure that the master branch will build and run.

This GitHub repository is designed to be an issue tracker, and a space for the public to submit patches through pull requests. It is synchronised with the main development repository that is available on the ICC's GitLab server which is available here.

Please feel free to submit issues to this repository, or even pull requests. We will try to deal with them as soon as possible, but as the core development team is quite small this could take some time.

Disclaimer

We would like to emphasise that SWIFT comes without any warranty of accuracy, correctness or efficiency. As mentioned in the license, the software comes as-is and the onus is on the user to get meaningful results. Whilst the authors will endeavour to answer questions related to using the code, we recommend users build and maintain their own copies. This documentation contains the most basic information to get started. Reading it and possibly also the source code is the best way to start running simulations.

The users are responsible to understand what the code is doing and for the results of their simulation runs.

Note also that the values of the parameters given in the examples are only indicative. We recommend users experiment by themselves and a campaign of experimentation with various values is highly encouraged. Each problem will likely require different values and the sensitivity to the details of the physical model is something left to the users to explore.

Acknowledgment & Citation

The SWIFT code was last described in this paper: https://ui.adsabs.harvard.edu/abs/2023arXiv230513380S. The core solver, the numerical methods as well as many extensions where described there. We ask users running SWIFT for their research to please cite this paper when they present their results.

In order to keep track of usage and measure the impact of the software, we kindly ask users publishing scientific results using SWIFT to add the following sentence to the acknowledgment section of their papers:

"The research in this paper made use of the SWIFT open-source simulation code (http://www.swiftsim.com, Schaller et al. 2018) version X.Y.Z."

with the version number set to the version used for the simulations and the reference pointing to the ASCL entry of the code: https://ascl.net/1805.020.

Contribution Guidelines

The SWIFT source code uses a variation of the 'Google' formatting style. The script 'format.sh' in the root directory applies the clang-format-18 tool with our style choices to all the SWIFT C source file. Please apply the formatting script to the files before submitting a pull request.

Please check that the test suite still runs with your changes applied before submitting a pull request and add relevant unit tests probing the correctness of new modules. An example of how to add a test to the suite can be found by considering the tests/testGreeting case.

Any contributions that fail any of the automated tests will not be accepted. Contributions that include tests of the proposed modules (or any current ones!) are highly encouraged.

Runtime parameters

 Welcome to the cosmological hydrodynamical code
    ______       _________________
   / ___/ |     / /  _/ ___/_  __/
   \__ \| | /| / // // /_   / /
  ___/ /| |/ |/ // // __/  / /
 /____/ |__/|__/___/_/    /_/
 SPH With Inter-dependent Fine-grained Tasking

 Version : 1.0.0
 Website: www.swiftsim.com
 Twitter: @SwiftSimulation

See INSTALL.swift for install instructions.

Usage: swift [options] [[--] param-file]
   or: swift [options] param-file
   or: swift_mpi [options] [[--] param-file]
   or: swift_mpi [options] param-file

Parameters:

    -h, --help                        show this help message and exit

  Simulation options:

    -b, --feedback                    Run with stars feedback.
    -c, --cosmology                   Run with cosmological time integration.
    --temperature                     Run with temperature calculation.
    -C, --cooling                     Run with cooling (also switches on --temperature).
    -D, --drift-all                   Always drift all particles even the ones
                                      far from active particles. This emulates
                                      Gadget-[23] and GIZMO's default behaviours.
    -F, --star-formation              Run with star formation.
    -g, --external-gravity            Run with an external gravitational potential.
    -G, --self-gravity                Run with self-gravity.
    -M, --multipole-reconstruction    Reconstruct the multipoles every time-step.
    -s, --hydro                       Run with hydrodynamics.
    -S, --stars                       Run with stars.
    -B, --black-holes                 Run with black holes.
    -k, --sinks                       Run with sink particles.
    -u, --fof                         Run Friends-of-Friends algorithm to
                                      perform black hole seeding.
    --lightcone                       Generate lightcone outputs.
    -x, --velociraptor                Run with structure finding.
    --line-of-sight                   Run with line-of-sight outputs.
    --limiter                         Run with time-step limiter.
    --sync                            Run with time-step synchronization
                                      of particles hit by feedback events.
    --csds                            Run with the Continuous Simulation Data
                                      Stream (CSDS).
    -R, --radiation                   Run with radiative transfer.
    --power                           Run with power spectrum outputs.

  Simulation meta-options:

    --quick-lyman-alpha               Run with all the options needed for the
                                      quick Lyman-alpha model. This is equivalent
                                      to --hydro --self-gravity --stars --star-formation
                                      --cooling.
    --eagle                           Run with all the options needed for the
                                      EAGLE model. This is equivalent to --hydro
                                      --limiter --sync --self-gravity --stars
                                      --star-formation --cooling --feedback
                                      --black-holes --fof.
    --gear                            Run with all the options needed for the
                                      GEAR model. This is equivalent to --hydro
                                      --limiter --sync --self-gravity --stars
                                      --star-formation --cooling --feedback.
    --agora                           Run with all the options needed for the
                                      GEAR model. This is equivalent to --hydro
                                      --limiter --sync --self-gravity --stars
                                      --star-formation --cooling --feedback.
                                      
  Control options:

    -a, --pin                         Pin runners using processor affinity.
    --nointerleave                    Do not interleave memory allocations across
                                      NUMA regions.
    -d, --dry-run                     Dry run. Read the parameter file, allocates
                                      memory but does not read the particles
                                      from ICs. Exits before the start of time
                                      integration. Checks the validity of
                                      parameters and IC files as well as memory
                                      limits.
    -e, --fpe                         Enable floating-point exceptions (debugging
                                      mode).
    -f, --cpu-frequency=<str>         Overwrite the CPU frequency (Hz) to be
                                      used for time measurements.
    -n, --steps=<int>                 Execute a fixed number of time steps.
                                      When unset use the time_end parameter
                                      to stop.
    -o, --output-params=<str>         Generate a parameter file with the options
                                      for selecting the output fields.
    -P, --param=<str>                 Set parameter value, overiding the value
                                      read from the parameter file. Can be used
                                      more than once {sec:par:value}.
    -r, --restart                     Continue using restart files.
    -t, --threads=<int>               The number of task threads to use on each
                                      MPI rank. Defaults to 1 if not specified.
    --pool-threads=<int>              The number of threads to use on each MPI
                                      rank for the threadpool operations.
                                      Defaults to the numbers of task threads
                                      if not specified.
    -T, --timers=<int>                Print timers every time-step.
    -v, --verbose=<int>               Run in verbose mode, in MPI mode 2 outputs
                                      from all ranks.
    -y, --task-dumps=<int>            Time-step frequency at which task graphs
                                      are dumped.
    --cell-dumps=<int>                Time-step frequency at which cell graphs
                                      are dumped.
    -Y, --threadpool-dumps=<int>      Time-step frequency at which threadpool
                                      tasks are dumped.
    --dump-tasks-threshold=<flt>      Fraction of the total step's time spent
                                      in a task to trigger a dump of the task plot
                                      on this step

See the file examples/parameter_example.yml for an example of parameter file.

soap's Issues

Precompute xray properties

It would be better to compute x-rays properties prior to running SOAP. This should be pretty easy to implement as SOAP already has the capacity for reading in extra properties for each particle, since that's how group_membership works.

SOAP is very slow on FLAMINGO high-res. high-z

Running the default script in the repo, the code times out after 72hours and seems to be very far from having completed.

SOAP fails on large runs due to halo catalogue arrays too large to MPI_GATHER()

Currently SOAP reads the input halo catalogue in parallel across all MPI ranks, but then gathers it on one rank to split it up into chunks. This fails on the L2800N10080 FLAMINGO run because it involves an allgather() where the count doesn't fit in an int.

The quick solution would be to implement our own gather which can handle large arrays using sends/recvs.

A more scalable solution would be to keep the halo catalogue distributed over all MPI ranks. When a chunk task is executed it would first have to request its assigned halos from whichever MPI rank they're stored on. Each rank would need to start a thread to respond to these requests. Since the chunk tasks are long running this thread could sleep most of the time and just check for messages periodically. The thread that hands out chunk tasks to nodes already works like this.

Projected velocity dispersions are very large

The following datasets have unreasonably large values in the FLAMINGO SOAP output:

/ProjectedAperture/??kpc/proj[xyz]/DarkMatterProjectedVelocityDispersion
/ProjectedAperture/??kpc/proj[xyz]/GasProjectedVelocityDispersion
/ProjectedAperture/??kpc/proj[xyz]/StellarProjectedVelocityDispersion

This seems to be due to a missing normalization factor that was added in 8c18a63 .

SOAP I/O optimization

SOAP seems to put a big load on the Lustre file system on Cosma so we aren't able to run many instances at the same time. With recent improvements in HDF5 we might be able to fix this.

SOAP splits the simulation volume into chunks and puts the chunks into a queue. When a compute node is idle it takes the next chunk from the queue and computes properties of the halos in the chunk. This continues until all chunks have been done. To read the particles for a chunk, the MPI ranks within the compute node all execute independent reads of the relevant parts of the input files. This is probably not very file system friendly.

We could instead have each node iterate through the files it needs to read and do a collective read on each one. This should use MPI-IO to funnel I/O through a few aggregator ranks, reducing the load on the file system while still having all ranks work on decompressing the data in HDF5.

Until recently this didn't work. HDF5 before 1.14 falls back to independent reads if you use any filters, and metadata operations always use independent I/O. But 1.14 can do collective reads and writes of filtered datasets and it also has an option to do metadata I/O using collective operations. If I use romio as the MPI-IO layer and enable lazy file opening (romio_no_indep_rw=true) then I can use collective HDF5 calls to read datasets while having only one rank per node actually touch the file system.

There's also a new function H5Dread_multi, which allows you to designate multiple datasets to read and which elements to read from them and then it reads everything in a single MPI-IO call. From some tests in C that seems to perform fairly well even with very few ranks accessing the file system. The down side is that H5Dread_multi doesn't exist in h5py and might be non-trivial to implement.

Rockstar support

We should try to read in Rockstar files and produce a catalog from there.

@kyleaoman mentioned he may have some old reading script. Maybe there is something there that can be used. Key bit is understanding the membership files and the subhalo <-> parent halo connection.

L1000N1800 group membership snapshot

The default snapshot number for the L1000N1800 group membership script is 78, however, there are only 77 snapshots for that run.
@jchelly

compute_halo_properties.py fails on snap8 due to "No space left on device"

Lots of runs failed before compute_halo_properties.py finished with the following message:

mca_fbtl_posix_pwritev: error in (p)write(v):No space left on device 
dynamic_gen2_write_all: fbtl_pwritev failed 

Depending on the point at which the run fails, SOAP produces halo_properties_XXXX.hdf5 with different file sizes. It is unclear which properties/datasets are incomplete in the final file. The compression does not fail on the next step.

The GAMA luminosities don't carry a unit

In the current version of SOAP, the GAMA luminosities defined for the stars don't carry a unit, while they could just carry erg/s. You can see this here:

SOAP/aperture_properties.py

Line 217 in 6100ac8

def stellar_luminosities(self):

Output ranking of subhalos by mass within their host halo

We would like to have the ranking of subhalos by mass within their host halo in the output. This is computed without using any particle data, unlike other SOAP halo properties.

The code does a collective read of the VR catalogue so we have the halos distributed across MPI ranks initially. We can compute the ranking by:

• parallel sort the subhalos by (host ID, descending total bound mass)
• identify sequences of subhalos in the same host locally on each MPI rank
• assign ranking by mass within each host locally on each MPI rank
• for the first host halo on each MPI rank we need to increment all subhalo ranks by the number of subhalos in that host on previous ranks in case the host's subhalos are split between MPI ranks

Then in halo_tasks.py we pass the new halo property through to the output.

Running SOAP not compatible with unyt 3.0.1

I have unyt 3.0.1 (the newest version)

When I git clone the master branch and run python3 -W error -m pytest aperture_properties.py

I get

[dc-brox1@login8a SOAP]$ python3 -W error -m pytest aperture_properties.py 
============================= test session starts ==============================
platform linux -- Python 3.9.1, pytest-7.4.0, pluggy-1.2.0
rootdir: /snap8/scratch/dp004/dc-brox1/SOAP_nov2023/SOAP
collected 1 item                                                               

aperture_properties.py F                                                 [100%]

=================================== FAILURES ===================================
___________________________ test_aperture_properties ___________________________

    def test_aperture_properties():
        """
        Unit test for the aperture property calculations.
    
        We generate 100 random "dummy" halos and feed them to
        ExclusiveSphereProperties::calculate() and
        InclusiveSphereProperties::calculate(). We check that the returned values
        are present, and have the right units, size and dtype
        """
    
        from dummy_halo_generator import DummyHaloGenerator
    
        # initialise the DummyHaloGenerator with a random seed
        dummy_halos = DummyHaloGenerator(3256)
        filter = RecentlyHeatedGasFilter(dummy_halos.get_cell_grid())
        stellar_age_calculator = StellarAgeCalculator(dummy_halos.get_cell_grid())
        cat_filter = CategoryFilter()
    
        pc_exclusive = ExclusiveSphereProperties(
            dummy_halos.get_cell_grid(), 50.0, filter, stellar_age_calculator, cat_filter
        )
        pc_inclusive = InclusiveSphereProperties(
            dummy_halos.get_cell_grid(), 50.0, filter, stellar_age_calculator, cat_filter
        )
    
        # generate 100 random halos
        for i in range(100):
            input_halo, data, _, _, _, particle_numbers = dummy_halos.get_random_halo(
                [1, 10, 100, 1000, 10000]
            )
            halo_result_template = {
                f"FOFSubhaloProperties/{PropertyTable.full_property_list['Ngas'][0]}": (
                    unyt.unyt_array(
                        particle_numbers["PartType0"],
                        dtype=PropertyTable.full_property_list["Ngas"][2],
                        units="dimensionless",
                    ),
                    "Dummy Ngas for filter",
                ),
                f"FOFSubhaloProperties/{PropertyTable.full_property_list['Ndm'][0]}": (
                    unyt.unyt_array(
                        particle_numbers["PartType1"],
                        dtype=PropertyTable.full_property_list["Ndm"][2],
                        units="dimensionless",
                    ),
                    "Dummy Ndm for filter",
                ),
                f"FOFSubhaloProperties/{PropertyTable.full_property_list['Nstar'][0]}": (
                    unyt.unyt_array(
                        particle_numbers["PartType4"],
                        dtype=PropertyTable.full_property_list["Nstar"][2],
                        units="dimensionless",
                    ),
                    "Dummy Nstar for filter",
                ),
                f"FOFSubhaloProperties/{PropertyTable.full_property_list['Nbh'][0]}": (
                    unyt.unyt_array(
                        particle_numbers["PartType5"],
                        dtype=PropertyTable.full_property_list["Nbh"][2],
                        units="dimensionless",
                    ),
                    "Dummy Nbh for filter",
                ),
            }
    
            for pc_type, pc_calc in [
                ("ExclusiveSphere", pc_exclusive),
                ("InclusiveSphere", pc_inclusive),
            ]:
                input_data = {}
                for ptype in pc_calc.particle_properties:
                    if ptype in data:
                        input_data[ptype] = {}
                        for dset in pc_calc.particle_properties[ptype]:
                            input_data[ptype][dset] = data[ptype][dset]
                input_halo_copy = input_halo.copy()
                input_data_copy = input_data.copy()
                halo_result = dict(halo_result_template)
>               pc_calc.calculate(input_halo, 0.0 * unyt.kpc, input_data, halo_result)

aperture_properties.py:1143: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
aperture_properties.py:956: in calculate
    part_props = ApertureParticleData(
aperture_properties.py:47: in __init__
    self.compute_basics()
aperture_properties.py:73: in compute_basics
    self.mass = unyt.array.uconcatenate(mass)
/cosma/home/dp004/dc-brox1/.local/lib/python3.9/site-packages/unyt/array.py:2291: in uconcatenate
    warn_deprecated(
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

name = 'unyt.uconcatenate'

    def warn_deprecated(
        name,
        /,
        *,
        stacklevel: int = 3,
        replacement: Optional[str] = None,
        since_version: str,
    ) -> None:
        msg = (
            f"{name} is deprecated and will be removed in a future version\n"
            f"Instead, {replacement}\n"
            f"(deprecated since unyt v{since_version})"
        )
>       warnings.warn(
            msg,
            category=DeprecationWarning,
            stacklevel=stacklevel,
        )
E       DeprecationWarning: unyt.uconcatenate is deprecated and will be removed in a future version
E       Instead, use numpy.concatenate
E       (deprecated since unyt v3.0)

/cosma/home/dp004/dc-brox1/.local/lib/python3.9/site-packages/unyt/_deprecation.py:18: DeprecationWarning
=========================== short test summary info ============================
FAILED aperture_properties.py::test_aperture_properties - DeprecationWarning: unyt.uconcatenate is deprecated and will be removed in ...
============================== 1 failed in 2.42s ==========

Compute centre of mass velocity of stars in apertures

Currently we compute the centre of mass position and velocity of the particles in each aperture. The apertures are centred on the potential minimum output by VR and we compute the mass weighted mean position and velocity of the particles within the aperture. All particle types are included in this calculation.

It would be useful to have the velocity of just the star particles too. We would compute the mass weighted mean velocity (and maybe also position?) of just the stars in each aperture and output that as a new halo property.

Possible issue with the compression of the compton-Y

From Scott:

Hi folks, is the L2800N5040 SOAP catalogue for snapshot 78 on COSMA8 provisional? I’m reading in Y parameters with infinite values. Currently extracting subset of key properties (centre, mass etc) for group/cluster haloes which are <1% of the total sample in this catalogue at z=0 (but still 2 million objects!)

@bwvdnbro's reaction:

That sounds like something went wrong with the compression. Unfortunately I (just) deleted the uncompressed version of that catalogue.

MassFractionSatellite calculation includes particles not in satellites

When calculating MassFractionSatellite we sum up the masses of all particles within the SO radius that are bound to a subhalo other than the central subhalo. However, we do not currently check if the subhalos which they are bound to are satellites of the halo. If you have two halos close to each other, then particles bound to halo 1 might be within the SO radius of halo 2. This leads to nonzero MassFractionSatellite values even if there are no satellites.

We want to implement a change so that MassFractionSatellite only considers particles bound to satellite subhalos. We will also include a new property, MassFractionExternal, which will capture the fraction of particles within the SO radius that are bound to other halos.

SOAP optimization: remove redundant particle sorting

SOAP could be sped up by feeding pre-sorted particles to the property calculations.

Currently the function process_single_halo() in halo_tasks.py finds all of the particles in some radius around the halo centre and feeds them to the property calculations from SO_properties.py, aperture_properties.py and subhalo_properties.py. The radius used is the largest radius needed by any of the property calculations. All property calculations receive the same set of particles, and all of the property calculations start by sorting the particles by radius.

There are two opportunities to speed up the code here:

• The sorting by radius could be done once before any of the property calculations are called
• Each property calculation could be passed just the particles it needs. Smaller apertures would then process many fewer particles.

To do this we would probably need to factor out the *ParticleData classes into a single class to be set up before calling any of the property calculations.

Document the group membership files

Before SOAP can be run we generate files with the halo membership for each particle in the snapshot. This is so that we can read sub-regions of the snapshot with halo membership information without needing to read the full VR output for every region processed.

These files may be useful for other purposes, so they should be documented.

SOAP could be packaged in a more standard way

This was provoked by some problems with module search paths while trying to run the halo matching script, which uses parts of SOAP.

Currently SOAP is just a bunch of source files in a directory and to run it you need to be in that directory or else you get import failures. We could turn it into a python package so that it could be pip installed (including dependencies) and we would no longer be relying on the current working directory to ensure all the modules are found.

We could also modify it to run with something like "python -m SOAP" instead of requiring the path to the main program source file.

Allow SOAP restarts

If SOAP crashes or times out we need to rerun it fully. We write out each chunks to a temporary file, so we could check if a chunk has already been calculated and skip if so

SOAP should write the bounding box of member particles

It would be useful to record the minimum and maximum position along each axis of particles in each halo - this can then be used directly to construct the spatial mask for a group. As it stands you'd have to take the centre, and try to guess what box around the centre is "big enough" to enclose all the particles in the list of member particles. No need to add bounding boxes for different radii, since those are trivially obtained from the centre and the radius - just the one for the unknown boundary of the member particles.

The only gotcha I can think of is the periodic boundary. I'd suggest that something close to the x boundary would have a bounding box like: [[0.9, 0.1], [0.4, 0.6], [0.4, 0.6]] (that's xmin, xmax, ymin, ymax, zmin, zmax). swiftsimio already knows to wrap the box when it sees xmin > xmax, so it can be used directly there, and in other contexts I think that the meaning is pretty intuitive (and if not easy to document).

compute_halo_properties.py fails for DMO runs.

compute_halo_properties.py fails for the latest subhalo ranking branch for DMO runs with the following message:

Traceback (most recent call last):
  File "/snap8/scratch/dp004/fkgm22/SOAP//compute_halo_properties.py", line 497, in <module>
    compute_halo_properties()
  File "/snap8/scratch/dp004/fkgm22/SOAP//compute_halo_properties.py", line 433, in compute_halo_properties
    metadata = task_queue.execute_tasks(
  File "/snap8/scratch/dp004/fkgm22/SOAP/task_queue.py", line 195, in execute_tasks
    result.append(task(*args))
  File "/snap8/scratch/dp004/fkgm22/SOAP/chunk_tasks.py", line 265, in __call__
    total_time, task_time, nr_left, nr_done = process_halos(comm, cellgrid.snap_unit_registry, data, mesh,
  File "/snap8/scratch/dp004/fkgm22/SOAP/halo_tasks.py", line 300, in process_halos
    halo_result = process_single_halo(
  File "/snap8/scratch/dp004/fkgm22/SOAP/halo_tasks.py", line 116, in process_single_halo
    idx_he, idx_T, idx_n, t_z, d_z, t_T, d_T, t_nH, d_nH, t_He, d_He, abundance_to_solar, joint_mask, volumes, data_n = xray_calc.find_indices(particle_data[ptype]['Densities'], particle_data[ptype]['Temperatures'], particle_data[ptype]['SmoothedElementMassFractions'], particle_data[ptype]['Masses'], fill_value = 0)
KeyError: 'PartType0'

The submission script can be found here:

/snap8/scratch/dp004/fkgm22/SOAP/scripts/FLAMINGO/L1000N1800

@Joeybraspenning, it might be related to the addition of the xray_calculator

Move to numpy 2.0

I've fixed the numpy version in the requirements file to < 2 for the time being

Bug in SO_properties.py

When I try to run python3 -W error -m pytest SO_properties.py

the code crashes with the error

=================================== FAILURES ===================================
______________________________ test_SO_properties ______________________________

    def test_SO_properties():
        from dummy_halo_generator import DummyHaloGenerator
    
        dummy_halos = DummyHaloGenerator(4251)
        filter = RecentlyHeatedGasFilter(dummy_halos.get_cell_grid())
        cat_filter = CategoryFilter()
    
        property_calculator_50kpc = SOProperties(
            dummy_halos.get_cell_grid(), filter, cat_filter, 50.0, "physical"
        )
        property_calculator_2500mean = SOProperties(
            dummy_halos.get_cell_grid(), filter, cat_filter, 2500.0, "mean"
        )
        property_calculator_2500crit = SOProperties(
            dummy_halos.get_cell_grid(), filter, cat_filter, 2500.0, "crit"
        )
        property_calculator_BN98 = SOProperties(
            dummy_halos.get_cell_grid(), filter, cat_filter, 0.0, "BN98"
        )
        property_calculator_5x2500mean = RadiusMultipleSOProperties(
            dummy_halos.get_cell_grid(), filter, cat_filter, 2500.0, 5.0, "mean"
        )
    
        for i in range(100):
            (
                input_halo,
                data,
                rmax,
                Mtot,
                Npart,
                particle_numbers,
            ) = dummy_halos.get_random_halo([2, 10, 100, 1000, 10000], has_neutrinos=True)
            halo_result_template = {
                f"FOFSubhaloProperties/{PropertyTable.full_property_list['Ngas'][0]}": (
                    unyt.unyt_array(
                        particle_numbers["PartType0"],
                        dtype=PropertyTable.full_property_list["Ngas"][2],
                        units="dimensionless",
                    ),
                    "Dummy Ngas for filter",
                ),
                f"FOFSubhaloProperties/{PropertyTable.full_property_list['Ndm'][0]}": (
                    unyt.unyt_array(
                        particle_numbers["PartType1"],
                        dtype=PropertyTable.full_property_list["Ndm"][2],
                        units="dimensionless",
                    ),
                    "Dummy Ndm for filter",
                ),
                f"FOFSubhaloProperties/{PropertyTable.full_property_list['Nstar'][0]}": (
                    unyt.unyt_array(
                        particle_numbers["PartType4"],
                        dtype=PropertyTable.full_property_list["Nstar"][2],
                        units="dimensionless",
                    ),
                    "Dummy Nstar for filter",
                ),
                f"FOFSubhaloProperties/{PropertyTable.full_property_list['Nbh'][0]}": (
                    unyt.unyt_array(
                        particle_numbers["PartType5"],
                        dtype=PropertyTable.full_property_list["Nbh"][2],
                        units="dimensionless",
                    ),
                    "Dummy Nbh for filter",
                ),
            }
            rho_ref = Mtot / (4.0 / 3.0 * np.pi * rmax ** 3)
    
            # force the SO radius to be outside the search sphere and check that
            # we get a ReadRadiusTooSmallError
            property_calculator_2500mean.reference_density = 0.01 * rho_ref
            property_calculator_2500crit.reference_density = 0.01 * rho_ref
            property_calculator_BN98.reference_density = 0.01 * rho_ref
            for prop_calc in [
                property_calculator_2500mean,
                property_calculator_2500crit,
                property_calculator_BN98,
            ]:
                fail = False
                try:
                    halo_result = dict(halo_result_template)
                    prop_calc.calculate(input_halo, rmax, data, halo_result)
                except ReadRadiusTooSmallError:
                    fail = True
                # 1 particle halos don't fail, since we always assume that the first
                # particle is at the centre of potential (which means we exclude it
                # in the SO calculation)
                # non-centrals don't fail, since we do not calculate any SO
                # properties and simply return zeros in this case
                assert (Npart == 1) or input_halo["Structuretype"] != 10 or fail
    
            # force the radius multiple to trip over not having computed the
            # required radius
            fail = False
            try:
                halo_result = dict(halo_result_template)
                property_calculator_5x2500mean.calculate(
                    input_halo, rmax, data, halo_result
                )
            except RuntimeError:
                fail = True
            assert fail
    
            # force the radius multiple to trip over the search radius
            fail = False
            try:
                halo_result = dict(halo_result_template)
                halo_result.update(
                    {
                        f"SO/2500_mean/{property_calculator_5x2500mean.radius_name}": (
                            0.1 * rmax,
                            "Dummy value.",
                        )
                    }
                )
                property_calculator_5x2500mean.calculate(
                    input_halo, 0.2 * rmax, data, halo_result
                )
            except ReadRadiusTooSmallError:
                fail = True
            assert fail
    
            # force the SO radius to be within the search sphere
            property_calculator_2500mean.reference_density = 2.0 * rho_ref
            property_calculator_2500crit.reference_density = 2.0 * rho_ref
            property_calculator_BN98.reference_density = 2.0 * rho_ref
    
            for SO_name, prop_calc in [
                ("50_kpc", property_calculator_50kpc),
                ("2500_mean", property_calculator_2500mean),
                ("2500_crit", property_calculator_2500crit),
                ("BN98", property_calculator_BN98),
                ("5xR_2500_mean", property_calculator_5x2500mean),
            ]:
                halo_result = dict(halo_result_template)
                # make sure the radius multiple is found this time
                if SO_name == "5xR_2500_mean":
                    halo_result[
                        f"SO/2500_mean/{property_calculator_5x2500mean.radius_name}"
                    ] = (0.1 * rmax, "Dummy value to force correct behaviour")
                input_data = {}
                for ptype in prop_calc.particle_properties:
                    if ptype in data:
                        input_data[ptype] = {}
                        for dset in prop_calc.particle_properties[ptype]:
                            input_data[ptype][dset] = data[ptype][dset]
                input_halo_copy = input_halo.copy()
                input_data_copy = input_data.copy()
>               prop_calc.calculate(input_halo, rmax, input_data, halo_result)

SO_properties.py:1916: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
SO_properties.py:1636: in calculate
    val = getattr(part_props, name)
lazy_properties.py:18: in _lazy_property
    setattr(self, attr_name, fn(self))
SO_properties.py:853: in Xraylum_restframe
    return self.gas_xraylum_restframe.sum(axis=0)
lazy_properties.py:18: in _lazy_property
    setattr(self, attr_name, fn(self))
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = <SO_properties.SOParticleData object at 0x2b3498df1e80>

    @lazy_property
    def gas_xraylum_restframe(self):
        if self.Ngas == 0:
            return None
>       return self.data["PartType0"]["XrayLuminositiesRestframe"][self.gas_selection]
E       KeyError: 'XrayLuminositiesRestframe'

SO_properties.py:841: KeyError
=========================== short test summary info ============================
FAILED SO_properties.py::test_SO_properties - KeyError: 'XrayLuminositiesRestframe'
============================== 1 failed in 2.57s ===

`unyt` breaking changes

I have been wrestling with SOAP issues, and realised that the most recent unyt version introduces breaking changes, such as breaking zero comparison, as in e.g.

SOAP/mask_cells.py

Line 30 in 0aec352

if done[halo_nr] == 0:

in unyt v2.9.2 this comparison yields True for a 0 dimensionless entry, but False in v2.9.4

This can be seen using unyt alone, first with unyt 2.9.4:

>>> unyt.__version__
'v2.9.4'
>>> a = unyt.unyt_array(np.zeros(5, dtype=np.int8), dtype=np.int8)
>>> for i in a:
...     print(i == 0)
... 
False
False
False
False
False
>>> 

Then with unyt 2.9.2 (latest previous version tested):

>>> unyt.__version__
'v2.9.2'
>>> a = unyt.unyt_array(np.zeros(5, dtype=np.int8), dtype=np.int8)
>>> for i in a:
...     print(i == 0)
... 
True
True
True
True
True
>>> 

comparing i.value == 0 yields Trues in both cases, but my short term suggestion is to modify requirements.txt to limit unyt<=2.9.2 as I don't know how pervasive the issue is (and know from testing this zero comparison issue comes up elsewhere).

Documentation table for SO

In the big table in the documentation the SO column always shows an 'X' for every single property and never a checkmark.
The checkmark should be set in the line of code I link below, but for some reason it doesn't work. I tried to figure out why not, but can't seem to find anything wrong.

https://github.com/SWIFTSIM/SOAP/blob/0a9400f67c8adb6c5dfdd65b65ed9b0517f49975/property_table.py#L1564C12-L1564C12

Add halo concentration following Wang+23

https://ui.adsabs.harvard.edu/abs/2023arXiv231000200W/abstract equations 3 and 5

Bug in io_test.py

When running mpirun -np 8 io_test.py the code crashes with error

Traceback (most recent call last):
  File "/cosma/home/dp004/dc-mcgi1/pull_requests/SOAP_32/io_test.py", line 82, in <module>
    io_test()
  File "/cosma/home/dp004/dc-mcgi1/pull_requests/SOAP_32/io_test.py", line 57, in io_test
    mesh = shared_mesh.SharedMesh(
  File "/cosma/home/dp004/dc-mcgi1/pull_requests/SOAP_32/shared_mesh.py", line 48, in __init__
    comm.Allreduce(pos_min_local, self.pos_min, op=MPI.MIN)
  File "mpi4py/MPI/Comm.pyx", line 876, in mpi4py.MPI.Comm.Allreduce
Traceback (most recent call last):
  File "/cosma/home/dp004/dc-mcgi1/pull_requests/SOAP_32/io_test.py", line 82, in <module>
  File "mpi4py/MPI/msgbuffer.pxi", line 748, in mpi4py.MPI._p_msg_cco.for_allreduce
    io_test()
  File "/cosma/home/dp004/dc-mcgi1/pull_requests/SOAP_32/io_test.py", line 57, in io_test
  File "mpi4py/MPI/msgbuffer.pxi", line 701, in mpi4py.MPI._p_msg_cco.for_cro_recv
    mesh = shared_mesh.SharedMesh(
  File "/cosma/home/dp004/dc-mcgi1/pull_requests/SOAP_32/shared_mesh.py", line 48, in __init__
  File "mpi4py/MPI/msgbuffer.pxi", line 203, in mpi4py.MPI.message_simple
    comm.Allreduce(pos_min_local, self.pos_min, op=MPI.MIN)
  File "mpi4py/MPI/Comm.pyx", line 876, in mpi4py.MPI.Comm.Allreduce
  File "mpi4py/MPI/msgbuffer.pxi", line 145, in mpi4py.MPI.message_basic
  File "mpi4py/MPI/msgbuffer.pxi", line 748, in mpi4py.MPI._p_msg_cco.for_allreduce
KeyError: '<d'
  File "mpi4py/MPI/msgbuffer.pxi", line 701, in mpi4py.MPI._p_msg_cco.for_cro_recv
  File "mpi4py/MPI/msgbuffer.pxi", line 203, in mpi4py.MPI.message_simple
  File "mpi4py/MPI/msgbuffer.pxi", line 145, in mpi4py.MPI.message_basic
KeyError: '<d'

We should also add a script that runs all tests

Use a parameter file instead of command line arguments

The scripts to run SOAP have become complicated and difficult to use, especially with the changes in #35 which allow use of different halo finders. It would probably be better to define the various paths and flags in a yaml parameter file.

Bug in property_table.py

When I try to run a control test of the SOAP code python3 -W error -m pytest property_table.py the test fails.

I get the error

______________________ ERROR collecting property_table.py ______________________
/cosma/home/dp004/dc-brox1/.local/lib/python3.9/site-packages/_pytest/python.py:617: in _importtestmodule
    mod = import_path(self.path, mode=importmode, root=self.config.rootpath)
/cosma/home/dp004/dc-brox1/.local/lib/python3.9/site-packages/_pytest/pathlib.py:565: in import_path
    importlib.import_module(module_name)
/cosma/local/Python/3.9.1-cosma8/lib/python3.9/importlib/__init__.py:127: in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
<frozen importlib._bootstrap>:1030: in _gcd_import
    ???
<frozen importlib._bootstrap>:1007: in _find_and_load
    ???
<frozen importlib._bootstrap>:986: in _find_and_load_unlocked
    ???
<frozen importlib._bootstrap>:680: in _load_unlocked
    ???
/cosma/home/dp004/dc-brox1/.local/lib/python3.9/site-packages/_pytest/assertion/rewrite.py:169: in exec_module
    source_stat, co = _rewrite_test(fn, self.config)
/cosma/home/dp004/dc-brox1/.local/lib/python3.9/site-packages/_pytest/assertion/rewrite.py:351: in _rewrite_test
    tree = ast.parse(source, filename=strfn)
/cosma/local/Python/3.9.1-cosma8/lib/python3.9/ast.py:50: in parse
    return compile(source, filename, mode, flags,
E     File "/snap8/scratch/dp004/dc-brox1/SOAP_nov2023/SOAP/property_table.py", line 35
E       output.append("\-")
E                     ^
E   SyntaxError: invalid escape sequence \-
=========================== short test summary info ============================
ERROR property_table.py

I did not change this function. The code does seem to run properly when I change "\-" to r"\-" in line 35

Filter aperture properties based on virial radius

We don't need to calculate large apertures for small halos. We could estimate the virial radius of each halo based on its bound mass and redshift. If an aperture is greater than some multiple N of the virial radius we would skip calculating it. The value of N would be set in the parameter file.

Colibre parameter file out of date?

I think the parameter file for Colibre might be out of date: parameter_files/colibre_SOAP_params.yml.

The latest version of the pipeline tries to use quantities that seem to not be saved in the SOAP catalogues. For example:

'CatalogueGroup_ExclusiveSphere_50kpc' object has no attribute 'gasmassincolddensediffusemetals'

In the parameter file itself, this is explicitly set to not be recorded:
GasMassInColdDenseDiffuseMetals: false

I suspect there are other fields that are inconsistent.

Remove FOFSubhaloProperties code

Now that we have moved away from VR we should remove the FOFSubhaloProperties calculation code. There is no likely use case for 6D FOF properties, and the name FOF was assumed by many people to be a 3D FOF. If we ever need to calculate 6D VR FOF quantities then we can just write a read routine that takes in all the bound+unbound particles

Ranks pausing when waiting for distributed input catalogue

See #48