doc/onboardingGuide/source/runtime_options.rst

+.. Runtime Options
+   Josh Borrow, 5th April 2018
+
+Runtime Options and Parameter Files
+===================================
+
+SWIFT requires a number of runtime options to run and get any sensible output.
+For instance, just running the ``swift`` binary will not use any SPH or gravity;
+the particles will just sit still!
+
+A list of  command line options can be found by running the compiled binary with
+the ``-h`` or ``--help`` flag:
+
+.. code-block:: bash
+
+   ./swift --help
+
+
+You will also need to specify a number of runtime parameters that are dependent 
+on your compile-time configuration in a parameter file. A list of all of these 
+parameters can be found in ``examples/parameter_example.yml``, and you can check 
+out examples in the ``examples/`` directory.

doc/onboardingGuide/source/submission_script.rst

+.. Submission Script
+
+Submission Script
+=================
+
+
+Below is an example submission script for the SLURM batch system. This runs 
+SWIFT with MPI, thread pinning, hydrodynamics, and self-gravity.
+
+.. code-block:: bash
+
+    #SBATCH --partition=<queue>
+    #SBATCH --account-name=<groupName>
+    #SBATCH --job-name=<jobName>
+    #SBATCH --nodes=<nNodes>
+    #SBATCH --ntasks-per-node=<nMPIRank>
+    #SBATCH --cpus-per-task=<nThreadsPerMPIRank>
+    #SBATCH --time=<hh>:<mm>:<ss>
+
+    srun -n $SLURM_NTASKS ./swift_mpi \ 
+        --threads=$SLURM_CPUS_PER_TASK --pin \
+        --hydro --self-gravity parameter_file.yml
+

doc/onboardingGuide/source/useful_configuration_flags.rst

+.. Configuration Options
+   Josh Borrow, 5th April 2018
+
+Useful Configuration Flags
+==========================
+
+A description of the available options for all flags including the examples
+below can be found by using
+``./configure  --help``.
+
+``--with-hydro=sphenix``
+~~~~~~~~~~~~~~~~~~~~~~~~
+There are several hydrodynamical schemes available in SWIFT. You can choose
+between them at compile-time with this option.
+
+``--with-riemann-solver=none``
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Some hydrodynamical schemes, for example GIZMO, require a Riemann solver.
+
+``--with-kernel=cubic-spline``
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Several kernels are made available for use with the hydrodynamical schemes.
+Choose between them with this compile-time flag.
+
+``--with-hydro-dimension=3``
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Run problems in 1, 2, and 3 (default) dimensions.
+
+``--with-equation-of-state=ideal-gas``
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Several equations of state are made available with this flag.
+
+``--with-cooling=none``
+~~~~~~~~~~~~~~~~~~~~~~~
+Several cooling implementations (including GRACKLE) are available.
+
+``--with-ext-potential=none``
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Many external potentials are available for use with SWIFT. 
+
+.. You can choose
+   between them at compile time. Some examples include a central potential, a
+   softened central potential, and a sinusoidal potential. You will need to
+   configure, for example, the mass in your parameter file at runtime.
+
+

doc/onboardingGuide/source/what_swift_can_do.rst

+What SWIFT can do
+=================
+
+SWIFT can solve a variety of problems aimed at cosmological and
+astrophysical applications. SWIFT's features include:
+
++ Hydrodynamics, using a variety of particle methods
++ Planetary science, with e.g. multiple equations of state
++ Dark Matter
++ Neutrinos
++ Gravity: self-gravity and external potentials
++ Cosmology 
++ Radiative cooling
++ Radiative transfer
++ On-the-fly analysis: halo finding (FOF), power spectra
++ And more!
+
+To enable and use these features, SWIFT needs to be compiled accordingly
+and corresponding flags need to be passed at runtime. Please consult the 
+instructions provided in the `documentation <https://swiftsim.com/docs>`_ 
+for full details.

examples/AGORA/AgoraDisk/README

+To run this example, SWIFT must be configured with the following options:
+
+./configure --with-subgrid=AGORA
+
+or similarily 
+
+./configure --with-feedback=AGORA --with-chemistry=AGORA  --with-cooling=grackle_0  --with-pressure-floor=GEAR --with-stars=GEAR --with-star-formation=GEAR 
+
+
+An analysis script for this test case can be found in the AGORA project repository:
+https://bitbucket.org/mornkr/agora-analysis-script/
+
+A modified version of the script able to read in SWIFT data is available upon request.
+
+(Note to the SWIFT authors: The modified script is stored with the other reference solutions.)

examples/AGORA/AgoraDisk/agora_disk.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e21 # kpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
+
+Scheduler:
+  max_top_level_cells: 16
+  cell_extra_sparts:  5000       # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
+  cell_extra_gparts:  5000       # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
+  cell_max_size:             8000   # (Optional) Maximal number of interactions per task if we force the split (this is the default value).
+  cell_sub_size_pair_hydro:  2560 # (Optional) Maximal number of hydro-hydro interactions per sub-pair hydro/star task (this is the default value).
+  cell_sub_size_self_hydro:  3200     # (Optional) Maximal number of hydro-hydro interactions per sub-self hydro/star task (this is the default value).
+  cell_sub_size_pair_stars:  2560 # (Optional) Maximal number of hydro-star interactions per sub-pair hydro/star task (this is the default value).
+  cell_sub_size_self_stars:  3200     # (Optional) Maximal number of hydro-star interactions per sub-self hydro/star task (this is the default value).
+  cell_sub_size_pair_grav:   2560 # (Optional) Maximal number of interactions per sub-pair gravity task  (this is the default value).
+  cell_sub_size_self_grav:   3200     # (Optional) Maximal number of interactions per sub-self gravity task  (this is the default value).
+  cell_split_size:           100       # (Optional) Maximal number of particles per cell (this is the default value).
+  
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   0.25  # The end time of the simulation (in internal units).
+  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
+  dt_max:     0.1  # The maximal time-step size of the simulation (in internal units).
+  max_dt_RMS_factor: 0.25  # (Optional) Dimensionless factor for the maximal displacement allowed based on the RMS velocities.
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            snapshot # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
+  compression:         4
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-3 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                       0.05    # Constant dimensionless multiplier for time integration.
+  MAC:                       geometric
+  theta_cr:                  0.7     
+  use_tree_below_softening:  1
+  comoving_DM_softening:     0.08 # Comoving softening length (in internal units).
+  max_physical_DM_softening: 0.08    # Physical softening length (in internal units).
+  comoving_baryon_softening:     0.08 # Comoving softening length (in internal units).
+  max_physical_baryon_softening: 0.08    # Physical softening length (in internal units).
+  mesh_side_length:       32        # Number of cells along each axis for the periodic gravity mesh.
+  
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+  h_min_ratio:           0.1095       # (Optional) Minimal allowed smoothing length in units of the softening. Defaults to 0 if unspecified.
+  h_max:                     10.      # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./agora_disk.hdf5     # The file to read
+  periodic:   0                     # Non-periodic BCs
+  cleanup_h_factors: 1              # Remove the h-factors inherited from Gadget
+  shift:    [674.1175, 674.1175, 674.1175]   # Centre the box
+
+
+# Cooling with Grackle 2.0
+GrackleCooling:
+  cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+  with_UV_background: 1 # Enable or not the UV background
+  redshift: 0 # Redshift to use (-1 means time based redshift)
+  with_metal_cooling: 1 # Enable or not the metal cooling
+  provide_volumetric_heating_rates: 0 # User provide volumetric heating rates
+  provide_specific_heating_rates: 0 # User provide specific heating rates
+  self_shielding_method: -1 # Grackle (<= 3) or Gear self shielding method
+  self_shielding_threshold_atom_per_cm3: 1e10  # Required only with GEAR's self shielding. Density threshold of the self shielding
+  max_steps: 1000
+  convergence_limit: 1e-2
+  thermal_time_myr: 0
+  maximal_density_Hpcm3: -1   # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter.
+
+
+GEARStarFormation:
+  star_formation_mode: agora        # default or agora
+  star_formation_efficiency: 0.01   # star formation efficiency (c_*)
+  maximal_temperature_K:     1e10   # Upper limit to the temperature of a star forming particle
+  density_threshold_Hpcm3:   10     # Density threshold in Hydrogen atoms/cm3
+  n_stars_per_particle: 1
+  min_mass_frac: 0.5
+
+
+GEARPressureFloor:
+  jeans_factor: 10
+
+AGORAFeedback:
+  energy_in_erg_per_CCSN: 1e51
+  supernovae_efficiency: 1
+  supernovae_explosion_time_myr: 5 # In Myr
+  ccsne_per_solar_mass : 0.010989 # 1/91
+  ejected_mass_in_solar_mass_per_CCSN : 14.8
+  ejected_Fe_mass_in_solar_mass_per_CCSN : 2.63
+  ejected_metal_mass_in_solar_mass_per_CCSN : 2.63
+
+AGORAChemistry:
+  initial_metallicity: 1          # if less than 0, read the metallicity from the snapshot (1 means solar abundance)
+  solar_abundance_Metals: 0.02  
+  scale_initial_metallicity: 1    # scale the initial metallicity using solar_abundance_Metals ?
+
+
+
+
+
+Restarts:
+  delta_hours:   1.        # (Optional) decimal hours between dumps of restart files.
+

examples/AGORA/AgoraDisk/checkSolution.py

+import h5py
+import numpy as np
+
+expected_total_mass = 227413100.0
+
+f = h5py.File("snapshot_0025.hdf5", "r")
+
+# get the total stellar mass
+mass = f["StarsParticles/Masses"]
+total_mass = sum(mass) * 1e10
+
+error = np.fabs((expected_total_mass - total_mass) / expected_total_mass)
+
+if error < 0.01:
+    print(48 * "!")
+    print(r"Well done !")
+    print(r"The stellar mass is %g Msol," % (total_mass))
+    print(r"The expected stellar mass is %g Msol" % (expected_total_mass))
+    print(r"This represents an error of %d %% !" % (100 * error))
+    print(48 * "!")
+else:
+    print(48 * "!")
+    print(r"Too bad !")
+    print(r"The stellar mass is %g Msol," % (total_mass))
+    print(r"While the expected stellar mass is %g Msol" % (expected_total_mass))
+    print(r"This represents an error of %d %% !" % (100 * error))
+    print(r"This is beyond the requirements.")
+    print(48 * "!")

examples/AGORA/AgoraDisk/getGrackleCoolingTable.sh

+#!/bin/bash
+
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5

examples/AGORA/AgoraDisk/getIC.sh

+#!/bin/bash
+
+if [ "$#" -ne 1 ]; then
+    echo "You need to provide the resolution (e.g. ./getIC.sh LOW)."
+    echo "The possible options are LOW and MED."
+    exit
+fi
+
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/snap.$1.hdf5
+mv snap.$1.hdf5 agora_disk.hdf5

examples/AGORA/AgoraDisk/run.sh

+#!/bin/bash
+
+# This example is based on the AGORA disk article (DOI: 10.3847/1538-4357/833/2/202)
+
+# set the resolution (LOW or MED)
+sim=LOW
+
+
+# make run.sh fail if a subcommand fails
+set -e
+
+# Generate the initial conditions if they are not present.
+if [ ! -e agora_disk.hdf5 ]
+then
+    echo "Fetching initial glass file for the Sedov blast example..."
+    ./getIC.sh $sim
+fi
+
+# Get the Grackle cooling table
+if [ ! -e CloudyData_UVB=HM2012.h5 ]
+then
+    echo "Fetching the Cloudy tables required by Grackle..."
+    ./getGrackleCoolingTable.sh
+fi
+
+
+
+
+# Run SWIFT
+../../../swift --sync --limiter --cooling --hydro --self-gravity --star-formation --feedback --stars --threads=8 agora_disk.yml 2>&1 | tee output.log
+
+
+
+

examples/ChemistryTests/MetalAdvectionTestEAGLE/README

+Metal advection test for the EAGLE chemistry scheme
+
+In some schemes, like the meshless (gizmo) scheme, particles exchange 
+masses via fluxes. This example tests whether the metals are correctly
+advected with those mass fluxes (for the EAGLE chemistry scheme)
+
+To run this test, compile with:
+    --with-hydro-dimension=2 --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-chemistry=EAGLE
+
+Due to the nature of this test, not much mass will be exchanged when running in Lagrangian mode,
+and hence not much advection will happen.
+To be able to see the effect of the advection, the hydrodynamics must be run in Eulerian mode.
+E.g. for gizmo-mvf: uncomment `#define GIZMO_FIX_PARTICLES` in src/const.h.
+
+Expected output when running in Eulerian mode is that the profiles should maintain their shape,
+but edges will be blurred due to numerical diffusion.

examples/ChemistryTests/MetalAdvectionTestEAGLE/advect_metals.yml

+MetaData:
+  run_name: advect_metals_EAGLE
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   0.1   # The end time of the simulation (in internal units).
+  dt_min:     1e-9  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-1  # The maximal time-step size of the simulation (in internal units).
+
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output # Common part of the name of output files
+  time_first:          0.     # Time of the first output (in internal units)
+  delta_time:          0.1    # Time between snapshots
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  time_first:          0.
+  delta_time:          1e-1   # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./advect_metals.hdf5     # The file to read
+  periodic:   1                      # periodic ICs.
+

examples/ChemistryTests/MetalAdvectionTestEAGLE/getGlass.sh

+#!/bin/bash
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_64.hdf5

examples/ChemistryTests/MetalAdvectionTestEAGLE/makeIC.py

+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Yolan Uyttenhove (yolan.uyttenhove@ugent.be)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+import math
+
+# ---------------------------------------------------------------------
+# Test the diffusion/advection of metals by creating regions with/without
+# initial metallicities. Run with EAGLE chemistry.
+# ---------------------------------------------------------------------
+
+import numpy as np
+import h5py
+
+GAMMA = 5 / 3
+RHO = 1
+P = 1
+VELOCITY = 2.5
+ELEMENT_COUNT = 9
+
+outputfilename = "advect_metals.hdf5"
+
+
+def get_masks(boxsize, pos):
+    masks = dict()
+    x = boxsize[0]
+    y = boxsize[1]
+
+    right_half_mask = pos[:, 0] > x / 2
+    masks["TotalMetallicity"] = right_half_mask
+    masks["He"] = pos[:, 0] * 2 % x > x / 2
+    masks["C"] = right_half_mask & (pos[:, 0] < 0.75 * x)
+    masks["N"] = right_half_mask & (pos[:, 0] > 0.75 * x)
+    masks["O"] = right_half_mask & (pos[:, 1] > 0.5 * y)
+    masks["Ne"] = right_half_mask & (pos[:, 1] < 0.5 * y)
+    masks["Mg"] = right_half_mask & (pos[:, 0] < 0.75 * x) & (pos[:, 1] > 0.5 * y)
+    masks["Si"] = right_half_mask & (pos[:, 0] > 0.75 * x) & (pos[:, 1] > 0.5 * y)
+    masks["Fe"] = right_half_mask & (pos[:, 0] < 0.75 * x) & (pos[:, 1] < 0.5 * y)
+
+    return masks
+
+
+def get_element_abundances_metallicity(pos, boxsize):
+    elements = ["He", "C", "N", "O", "Ne", "Mg", "Si", "Fe"]
+    abundances = [0.2, 0.2, 0.2, 0.2, 0.2, 0.1, 0.1, 0.1]
+
+    masks = get_masks(boxsize, pos)
+    total_metallicity = np.where(masks["TotalMetallicity"], 0.7, 0.1)
+    element_abundances = [
+        np.where(masks[k], v, 0) for k, v in zip(elements, abundances)
+    ]
+
+    # Finally add H so that H + He + TotalMetallicity = 1
+    return (
+        np.stack(
+            [1 - element_abundances[0] - total_metallicity] + element_abundances, axis=1
+        ),
+        total_metallicity,
+    )
+
+
+if __name__ == "__main__":
+    glass = h5py.File("glassPlane_64.hdf5", "r")
+    parts = glass["PartType0"]
+    pos = parts["Coordinates"][:]
+    pos = np.concatenate([pos, pos + np.array([1, 0, 0])])
+    h = parts["SmoothingLength"][:]
+    h = np.concatenate([h, h])
+    glass.close()
+
+    # Set up metadata
+    boxsize = np.array([2.0, 1.0, 1.0])
+    n_part = len(h)
+    ids = np.arange(n_part) + 1
+
+    # Setup other particle quantities
+    rho = RHO * np.ones_like(h)
+    rho[pos[:, 1] < 0.5 * boxsize[1]] *= 0.5
+    masses = rho * np.prod(boxsize) / n_part
+    velocities = np.zeros((n_part, 3))
+    velocities[:, :] = 0.5 * math.sqrt(2) * VELOCITY * np.array([1.0, 1.0, 0.0])
+    internal_energy = P / (rho * (GAMMA - 1))
+
+    # Setup metallicities
+    element_abundances, metallicities = get_element_abundances_metallicity(pos, boxsize)
+
+    # Now open the file and write the data.
+    file = h5py.File(outputfilename, "w")
+
+    # Header
+    grp = file.create_group("/Header")
+    grp.attrs["BoxSize"] = boxsize
+    grp.attrs["NumPart_Total"] = [n_part, 0, 0, 0, 0, 0]
+    grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+    grp.attrs["NumPart_ThisFile"] = [n_part, 0, 0, 0, 0, 0]
+    grp.attrs["Time"] = 0.0
+    grp.attrs["NumFilesPerSnapshot"] = 1
+    grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    grp.attrs["Flag_Entropy_ICs"] = 0
+    grp.attrs["Dimension"] = 2
+
+    # Units
+    grp = file.create_group("/Units")
+    grp.attrs["Unit length in cgs (U_L)"] = 1.0
+    grp.attrs["Unit mass in cgs (U_M)"] = 1.0
+    grp.attrs["Unit time in cgs (U_t)"] = 1.0
+    grp.attrs["Unit current in cgs (U_I)"] = 1.0
+    grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
+
+    # Particle group
+    grp = file.create_group("/PartType0")
+    grp.create_dataset("Coordinates", data=pos, dtype="d")
+    grp.create_dataset("Velocities", data=velocities, dtype="f")
+    grp.create_dataset("Masses", data=masses, dtype="f")
+    grp.create_dataset("SmoothingLength", data=h, dtype="f")
+    grp.create_dataset("InternalEnergy", data=internal_energy, dtype="f")
+    grp.create_dataset("ParticleIDs", data=ids, dtype="L")
+    grp.create_dataset("Metallicity", data=metallicities, dtype="f")
+    grp.create_dataset("ElementAbundance", data=element_abundances, dtype="f")
+
+    file.close()
+
+    # close up, and we're done!
+    file.close()

examples/ChemistryTests/MetalAdvectionTestEAGLE/plotAdvectedMetals.py

+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Yolan Uyttenhove (yolan.uyttenhove@ugent.be)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+import math
+
+import matplotlib.pyplot as plt
+from matplotlib.cm import ScalarMappable
+from matplotlib.colors import Normalize
+import numpy as np
+import unyt
+
+import swiftsimio
+from swiftsimio.visualisation import project_gas
+from pathlib import Path
+
+from makeIC import VELOCITY
+
+
+def plot_single(ax_mass, ax_fraction, name, title, data, mass_map, kwargs_inner):
+    mass_weighted_map = project_gas(data, project=name, **kwargs_inner["projection"])
+    ax_mass.imshow(mass_weighted_map.in_cgs().value.T, **kwargs_inner["imshow_mass"])
+    ax_mass.set_title(title)
+    ax_fraction.imshow(
+        (mass_weighted_map / mass_map).T, **kwargs_inner["imshow_fraction"]
+    )
+    ax_fraction.set_title(title)
+    ax_mass.axis("off")
+    ax_fraction.axis("off")
+
+
+def plot_all(fname, savename):
+    data = swiftsimio.load(fname)
+
+    # Shift coordinates to starting position
+    velocity = 0.5 * math.sqrt(2) * VELOCITY * np.array([1, 1, 0]) * unyt.cm / unyt.s
+    data.gas.coordinates -= data.metadata.time * velocity
+
+    # Add mass weighted element mass fractions to gas dataset
+    masses = data.gas.masses
+    element_mass_fractions = data.gas.element_mass_fractions
+    data.gas.element_mass_H = element_mass_fractions.hydrogen * masses
+    data.gas.element_mass_He = element_mass_fractions.helium * masses
+    data.gas.element_mass_C = element_mass_fractions.carbon * masses
+    data.gas.element_mass_N = element_mass_fractions.nitrogen * masses
+    data.gas.element_mass_O = element_mass_fractions.oxygen * masses
+    data.gas.element_mass_Ne = element_mass_fractions.neon * masses
+    data.gas.element_mass_Mg = element_mass_fractions.magnesium * masses
+    data.gas.element_mass_Si = element_mass_fractions.silicon * masses
+    data.gas.element_mass_Fe = element_mass_fractions.iron * masses
+    data.gas.total_metal_mass = data.gas.metal_mass_fractions * masses
+
+    # Create necessary figures and axes
+    fig = plt.figure(layout="constrained", figsize=(16, 8))
+    fig.suptitle(
+        f"Profiles shifted to starting position after t={data.metadata.time:.2f}",
+        fontsize=14,
+    )
+    fig_ratios, fig_masses = fig.subfigures(2, 1)
+    fig_ratios.suptitle("Mass ratio of elements")
+    fig_masses.suptitle("Surface density in elements")
+    axes_ratios = fig_ratios.subplots(2, 5, sharex=True, sharey=True)
+    axes_masses = fig_masses.subplots(2, 5, sharex=True, sharey=True)
+
+    # parameters for swiftsimio projections
+    projection_kwargs = {
+        "region": np.array([0, 2, 0, 1, 0, 1]) * unyt.cm,
+        "resolution": 500,
+        "parallel": True,
+    }
+    # Parameters for imshow
+    norm_ratios = Normalize(0, 0.95)
+    norm_masses = Normalize(0, 0.9)
+    imshow_fraction_kwargs = dict(norm=norm_ratios, cmap="rainbow")
+    imshow_mass_kwargs = dict(norm=norm_masses, cmap="turbo")
+
+    # Plot the quantities
+    mass_map = project_gas(data, project="masses", **projection_kwargs)
+
+    # Parameters for plotting:
+    plotting_kwargs = dict(
+        data=data,
+        mass_map=mass_map,
+        kwargs_inner=dict(
+            projection=projection_kwargs,
+            imshow_mass=imshow_mass_kwargs,
+            imshow_fraction=imshow_fraction_kwargs,
+        ),
+    )
+
+    plot_single(
+        axes_masses[0][0],
+        axes_ratios[0][0],
+        "total_metal_mass",
+        "All metals",
+        **plotting_kwargs,
+    )
+    plot_single(
+        axes_masses[0][1],
+        axes_ratios[0][1],
+        "element_mass_H",
+        "Hydrogen",
+        **plotting_kwargs,
+    )
+    plot_single(
+        axes_masses[0][2],
+        axes_ratios[0][2],
+        "element_mass_He",
+        "Helium",
+        **plotting_kwargs,
+    )
+    plot_single(
+        axes_masses[0][3],
+        axes_ratios[0][3],
+        "element_mass_C",
+        "Carbon",
+        **plotting_kwargs,
+    )
+    plot_single(
+        axes_masses[0][4],
+        axes_ratios[0][4],
+        "element_mass_N",
+        "Nitrogen",
+        **plotting_kwargs,
+    )
+    plot_single(
+        axes_masses[1][0],
+        axes_ratios[1][0],
+        "element_mass_O",
+        "Oxygen",
+        **plotting_kwargs,
+    )
+    plot_single(
+        axes_masses[1][1],
+        axes_ratios[1][1],
+        "element_mass_Ne",
+        "Neon",
+        **plotting_kwargs,
+    )
+    plot_single(
+        axes_masses[1][2],
+        axes_ratios[1][2],
+        "element_mass_Mg",
+        "Magnesium",
+        **plotting_kwargs,
+    )
+    plot_single(
+        axes_masses[1][3],
+        axes_ratios[1][3],
+        "element_mass_Si",
+        "Silicon",
+        **plotting_kwargs,
+    )
+    plot_single(
+        axes_masses[1][4],
+        axes_ratios[1][4],
+        "element_mass_Fe",
+        "Iron",
+        **plotting_kwargs,
+    )
+
+    # Add Colorbars
+    cb_masses = fig_masses.colorbar(
+        ScalarMappable(**imshow_mass_kwargs), shrink=0.75, pad=0.01, ax=axes_masses
+    )
+    cb_ratios = fig_ratios.colorbar(
+        ScalarMappable(**imshow_fraction_kwargs), shrink=0.75, pad=0.01, ax=axes_ratios
+    )
+    cb_masses.ax.set_ylabel("Surface density (g/cm^2)", rotation=270, labelpad=15)
+    cb_ratios.ax.set_ylabel("Mass ratio", rotation=270, labelpad=15)
+
+    # Save output
+    fig.savefig(savename, dpi=300)
+
+
+if __name__ == "__main__":
+    cwd = Path(__file__).parent
+    print("Plotting metals for output_0000.hdf5...")
+    plot_all(cwd / "output_0000.hdf5", savename=cwd / "metal_advection_0000.png")
+    print("Plotting metals for output_0001.hdf5...")
+    plot_all(cwd / "output_0001.hdf5", savename=cwd / "metal_advection_0001.png")
+    print("Done!")

examples/ChemistryTests/MetalAdvectionTestEAGLE/run.sh

+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -e glassPlane_64.hdf5 ]
+then
+    echo "Fetching initial glass file ..."
+    ./getGlass.sh
+fi
+
+if [ ! -f 'advect_metals.hdf5' ]; then
+    echo "Generating ICs"
+    python3 makeIC.py
+fi
+
+# Run SWIFT (must be compiled with --with-chemistry=EAGLE)
+../../../swift \
+    --hydro --threads=4 advect_metals.yml 2>&1 | tee output.log
+
+python3 runSanityChecksAdvectedMetals.py
+python3 plotAdvectedMetals.py

examples/ChemistryTests/MetalAdvectionTestEAGLE/runSanityChecksAdvectedMetals.py

+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Yolan Uyttenhove (yolan.uyttenhove@ugent.be)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+from pathlib import Path
+
+import numpy as np
+import swiftsimio
+
+
+def test(name):
+    def test_decorator(test_func):
+        def test_runner(*args, **kwargs):
+            try:
+                test_func(*args, **kwargs)
+            except Exception as e:
+                print(f"\tTest {name} failed!")
+                raise e
+            else:
+                print(f"\tTest {name} \u2705")
+
+        return test_runner
+
+    return test_decorator
+
+
+def approx_equal(a, b, threshold=0.001):
+    return 2 * abs(a - b) / (abs(a) + abs(b)) < threshold
+
+
+@test("mass fractions sum to 1")
+def test_sum_mass_fractions(data):
+    metal_fractions = data.gas.metal_mass_fractions
+    h_fraction = data.gas.element_mass_fractions.hydrogen
+    he_fraction = data.gas.element_mass_fractions.helium
+    total = metal_fractions + he_fraction + h_fraction
+
+    assert np.all(approx_equal(total, 1))
+
+
+@test("total metal mass conservation")
+def test_total_metal_mass_conservation(data_start, data_end):
+    def metal_mass(data):
+        return data.gas.masses * data.gas.metal_mass_fractions
+
+    assert approx_equal(np.sum(metal_mass(data_start)), np.sum(metal_mass(data_end)))
+
+
+def element_mass(data, element_name):
+    return data.gas.masses * getattr(data.gas.element_mass_fractions, element_name)
+
+
+@test("hydrogen mass conservation")
+def test_h_mass_conservation(data_start, data_end):
+    assert approx_equal(
+        np.sum(element_mass(data_start, "hydrogen")),
+        np.sum(element_mass(data_end, "hydrogen")),
+    )
+
+
+@test("helium mass conservation")
+def test_he_mass_conservation(data_start, data_end):
+    assert approx_equal(
+        np.sum(element_mass(data_start, "helium")),
+        np.sum(element_mass(data_end, "helium")),
+    )
+
+
+@test("carbon mass conservation")
+def test_c_mass_conservation(data_start, data_end):
+    assert approx_equal(
+        np.sum(element_mass(data_start, "carbon")),
+        np.sum(element_mass(data_end, "carbon")),
+    )
+
+
+@test("nitrogen mass conservation")
+def test_n_mass_conservation(data_start, data_end):
+    assert approx_equal(
+        np.sum(element_mass(data_start, "nitrogen")),
+        np.sum(element_mass(data_end, "nitrogen")),
+    )
+
+
+@test("oxygen mass conservation")
+def test_o_mass_conservation(data_start, data_end):
+    assert approx_equal(
+        np.sum(element_mass(data_start, "oxygen")),
+        np.sum(element_mass(data_end, "oxygen")),
+    )
+
+
+@test("neon mass conservation")
+def test_ne_mass_conservation(data_start, data_end):
+    assert approx_equal(
+        np.sum(element_mass(data_start, "neon")), np.sum(element_mass(data_end, "neon"))
+    )
+
+
+@test("magnesium mass conservation")
+def test_mg_mass_conservation(data_start, data_end):
+    assert approx_equal(
+        np.sum(element_mass(data_start, "magnesium")),
+        np.sum(element_mass(data_end, "magnesium")),
+    )
+
+
+@test("silicon mass conservation")
+def test_si_mass_conservation(data_start, data_end):
+    assert approx_equal(
+        np.sum(element_mass(data_start, "silicon")),
+        np.sum(element_mass(data_end, "silicon")),
+    )
+
+
+@test("iron mass conservation")
+def test_fe_mass_conservation(data_start, data_end):
+    assert approx_equal(
+        np.sum(element_mass(data_start, "iron")), np.sum(element_mass(data_end, "iron"))
+    )
+
+
+if __name__ == "__main__":
+    print("Running sanity checks...")
+
+    cwd = Path(__file__).parent
+    start = swiftsimio.load(cwd / "output_0000.hdf5")
+    end = swiftsimio.load(cwd / "output_0001.hdf5")
+
+    test_sum_mass_fractions(end)
+    test_total_metal_mass_conservation(start, end)
+    test_h_mass_conservation(start, end)
+    test_he_mass_conservation(start, end)
+    test_c_mass_conservation(start, end)
+    test_n_mass_conservation(start, end)
+    test_o_mass_conservation(start, end)
+    test_ne_mass_conservation(start, end)
+    test_mg_mass_conservation(start, end)
+    test_si_mass_conservation(start, end)
+    test_fe_mass_conservation(start, end)
+
+    print("Done!")

examples/ChemistryTests/MetalAdvectionTestGEAR/README

+Metal advection test for the GEAR/AGORA chemistry schemes
+
+In some schemes, like the meshless (gizmo) scheme, particles exchange 
+masses via fluxes. This example tests whether the metals are correctly
+advected with those mass fluxes.
+
+To run this test with GEAR, compile with:
+    --with-hydro-dimension=2 --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-chemistry=GEAR_X
+with X the desired number of elements. The number of elements must be also be set in the makeIC.py script
+and must be the equal to one used in the compilation flag. The default value is 4.
+
+To run this test with AGORA, compile with:
+    --with-hydro-dimension=2 --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-chemistry=AGORA
+NOTE: The AGORA scheme only traces Fe mass and total metal mass, so the number of elements in makeIC.py must be set
+to the value of 2.
+
+Due to the nature of this test, not much mass will be exchanged when running in Lagrangian mode,
+and hence not much advection will happen.
+To be able to see the effect of the advection, the hydrodynamics must be run in Eulerian mode.
+E.g. for gizmo-mvf: uncomment `#define GIZMO_FIX_PARTICLES` in src/const.h.
+
+Expected output when running in Eulerian mode is that the profiles should maintain their shape,
+but edges will be blurred due to numerical diffusion.

examples/ChemistryTests/MetalAdvectionTestGEAR/advect_metals.yml

+MetaData:
+  run_name: advect_metals_GEAR
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   0.1   # The end time of the simulation (in internal units).
+  dt_min:     1e-9  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-1  # The maximal time-step size of the simulation (in internal units).
+
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output # Common part of the name of output files
+  time_first:          0.     # Time of the first output (in internal units)
+  delta_time:          0.1    # Time between snapshots
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  time_first:          0.
+  delta_time:          1e-1   # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./advect_metals.hdf5     # The file to read
+  periodic:   1                        # periodic ICs.
+
+GEARChemistry:
+  initial_metallicity: -1         # Negative values mean that the metallicity will be read from the ICs
+
+AGORAChemistry:
+  initial_metallicity: -1         # Negative values mean that the metallicity will be read from the ICs
+  scale_initial_metallicity: 1    # scale the initial metallicity using solar_abundance_Metals? (needed to prevent crash in writing of metadata)

examples/ChemistryTests/MetalAdvectionTestGEAR/getGlass.sh

+#!/bin/bash
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_64.hdf5