WIP: test for metal advection

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

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.5   # 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:          5e-2   # Time between snapshots
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  time_first:          0.
+  delta_time:          5e-2   # 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.9      # 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 http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_128.hdf5

examples/ChemistryTests/MetalAdvectionTestEAGLE/makeIC.py

+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Yolan Uyttenhove (yolan.uyttehove@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/>.
+#
+##############################################################################
+
+
+# ---------------------------------------------------------------------
+# 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 = 1
+BOX_SIZE = 1
+
+outputfilename = "advect_metals.hdf5"
+
+if __name__ == "__main__":
+
+    glass = h5py.File("glassPlane_128.hdf5", "r")
+    parts = glass["PartType0"]
+    pos = parts["Coordinates"][:]
+    h = parts["SmoothingLength"][:]
+    glass.close()
+
+    # Set up metadata
+    boxsize = np.array([1.0, 1.0, 0.0])
+    n_part = len(h)
+    ids = np.arange(n_part) + 1
+
+    # Setup other particle quantities
+    masses = RHO * BOX_SIZE ** 3 / n_part * np.ones(n_part)
+    velocities = np.zeros((n_part, 3))
+    velocities[:, 0] = VELOCITY
+    internal_energy = P / (RHO * (GAMMA - 1)) * np.ones(n_part)
+
+    # Setup metallicities
+    metallicities = np.zeros(n_part)
+    # Mask for middle square
+    mask = ((1 / 3 * BOX_SIZE < pos[:, 0]) & (pos[:, 0] < 2 / 3 * BOX_SIZE) &
+            (1 / 3 * BOX_SIZE < pos[:, 1]) & (pos[:, 1] < 2 / 3 * BOX_SIZE))
+    metallicities[mask] = 1
+
+    # Now open the file and write the data.
+    file = h5py.File(outputfilename, "w")
+
+    # Header
+    grp = file.create_group("/Header")
+    grp.attrs["BoxSize"] = [1.0, 1.0, 1.0]
+    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="L")
+
+    # TODO: add ElementAbundance arrays and IronMassFracFromSNIa (see EAGLE/chemistry_io.h)
+
+    file.close()
+
+    # close up, and we're done!
+    file.close()

examples/ChemistryTests/MetalAdvectionTestEAGLE/plotAdvectedMetals.py

+import matplotlib.pyplot as plt
+import sys
+
+import numpy as np
+import swiftsimio
+from swiftsimio.visualisation import project_gas
+import scipy.interpolate as si
+from pathlib import Path
+
+
+def project_nn(x, y, data, data_xy):
+    xv, yv = np.meshgrid(x, y)
+    return si.griddata(data_xy, data, (xv, yv), method="nearest")
+
+if __name__ == "__main__":
+    try:
+        n = sys.argv[1]
+    except IndexError:
+        n = 2
+
+    cwd = Path(__file__).parent
+
+    fname = cwd / f"output_{n:04}.hdf5"
+    # fname = cwd / f"IC.hdf5"
+
+    data = swiftsimio.load(fname)
+
+    # parameters for swiftsimio projections
+    projection_kwargs = {"resolution": 1024, "parallel": True}
+
+    mass_map = project_gas(data, project="masses", **projection_kwargs)
+    mass_weighted_metal_map = project_gas(data, project="metal_mass_fractions", **projection_kwargs)
+    metal_map = mass_weighted_metal_map / mass_map
+
+    res = 1000
+    plt.imshow(metal_map.T)
+    plt.axis("off")
+    plt.tight_layout()
+    plt.savefig("metal_advection.png", dpi=300)
+    plt.show()
\ No newline at end of file

examples/ChemistryTests/MetalAdvectionTestEAGLE/run.sh

+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -e glassPlane_128.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 plotAdvectedMetals.py

src/const.h

@@ -40,7 +40,7 @@
 
 /* Options to control the movement of particles for GIZMO_SPH. */
 /* This option disables particle movement */
-//#define GIZMO_FIX_PARTICLES
+#define GIZMO_FIX_PARTICLES
 /* Try to keep cells regular by adding a correction velocity. */
 //#define GIZMO_STEER_MOTION
 /* Use the total energy instead of the thermal energy as conserved variable. */