examples/Planetary/GreshoVortex_3D/gresho.yml → 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
@@ -6,38 +9,38 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
 
-Scheduler:
-  max_top_level_cells: 15
-
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   1.    # The end time of the simulation (in internal units).
-  dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
-  dt_max:     1e-2  # The maximal time-step size 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:            gresho # Common part of the name of output files
+  basename:            output # Common part of the name of output files
   time_first:          0.     # Time of the first output (in internal units)
-  delta_time:          1e-1   # Time difference between consecutive outputs (in internal units)
-  compression:         1
-  
+  delta_time:          0.1    # Time between snapshots
+
 # Parameters governing the conserved quantities statistics
 Statistics:
-  delta_time:          1e-2 # Time between statistics output
+  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.
-  viscosity_alpha:       0.8      # Override for the initial value of the artificial viscosity. In schemes that have a fixed AV, this remains as alpha throughout the run.
-  
+
 # Parameters related to the initial conditions
 InitialConditions:
-  file_name:  ./greshoVortex.hdf5     # The file to read
-  periodic:   1
+  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
 
-# Parameters related to the equation of state
-EoS:
-    planetary_use_idg_def:    1               # Default ideal gas, material ID 0
+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

examples/ChemistryTests/MetalAdvectionTestGEAR/makeIC.py

+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2024 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
+# Set ELEMENT_COUNT to 2 for AGORA runs, or for GEAR compile swift with `--with-chemistry=GEAR_X` with X equal to
+# ELEMENT_COUNT.
+ELEMENT_COUNT = 4
+
+outputfilename = "advect_metals.hdf5"
+
+
+def get_mask(boxsize, pos, element_idx):
+    x = boxsize[0]
+    right_half_mask = pos[:, 0] > x / 3
+
+    if element_idx == ELEMENT_COUNT - 1:
+        return right_half_mask
+    else:
+        periods = 2 * (element_idx + 2)
+        periods_mask = (pos[:, 0] * periods // x) % 2 == 0
+        return right_half_mask & periods_mask
+
+
+def get_abundance(element_idx):
+    if element_idx == ELEMENT_COUNT - 1:
+        return 0.2
+    else:
+        return 0.1 * 0.5 ** element_idx
+
+
+def get_element_abundances_metallicity(pos, boxsize):
+    element_abundances = []
+    for i in range(ELEMENT_COUNT):
+        element_abundances.append(
+            np.where(get_mask(boxsize, pos, i), get_abundance(i), 0)
+        )
+
+    return np.stack(element_abundances, axis=1)
+
+
+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
+    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("MetalMassFraction", data=metallicities, dtype="f")
+
+    file.close()
+
+    # close up, and we're done!
+    file.close()
+
+    if ELEMENT_COUNT == 2:
+        print(
+            f"Make sure swift was compiled with `--with-chemistry=GEAR_2` or `--with-chemistry=AGORA`"
+        )
+    else:
+        print(
+            f"Make sure swift was compiled with `--with-chemistry=GEAR_{ELEMENT_COUNT}`"
+        )

examples/ChemistryTests/MetalAdvectionTestGEAR/plotAdvectedMetals.py

+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2024 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 SymLogNorm, Normalize
+import numpy as np
+import unyt
+
+import swiftsimio
+from swiftsimio.visualisation import project_gas
+from pathlib import Path
+
+from makeIC import VELOCITY, ELEMENT_COUNT
+
+
+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).value.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.metal_mass_fractions
+    columns = getattr(element_mass_fractions, "named_columns", None)
+    for i in range(ELEMENT_COUNT):
+        if columns is None:
+            data.gas.__setattr__(
+                f"element_mass_sp{i}", element_mass_fractions[:, i] * masses
+            )
+        else:
+            data.gas.__setattr__(
+                f"element_mass_sp{i}",
+                getattr(element_mass_fractions, columns[i]) * masses,
+            )
+
+    # Create necessary figures and axes
+    fig = plt.figure(layout="constrained", figsize=(8, 2 * ELEMENT_COUNT + 1))
+    fig.suptitle(
+        f"Profiles shifted to starting position after t={data.metadata.time:.2f}",
+        fontsize=14,
+    )
+    fig_ratios, fig_masses = fig.subfigures(1, 2)
+    fig_ratios.suptitle("Mass ratio of elements")
+    fig_masses.suptitle("Surface density in elements")
+    axes_ratios = fig_ratios.subplots(ELEMENT_COUNT, 1, sharex=True, sharey=True)
+    axes_masses = fig_masses.subplots(ELEMENT_COUNT, 1, 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
+    if ELEMENT_COUNT > 5:
+        thresh = 10 ** math.floor(math.log10(0.5 ** (ELEMENT_COUNT - 2)))
+        norm_ratios = SymLogNorm(vmin=0, vmax=0.21, linthresh=thresh, base=10)
+        norm_masses = SymLogNorm(vmin=0, vmax=0.25, linthresh=thresh, base=10)
+    else:
+        norm_ratios = Normalize(vmin=0, vmax=0.21)
+        norm_masses = Normalize(vmin=0, vmax=0.25)
+    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,
+        ),
+    )
+
+    if columns is None:
+        columns = [f"Species {i + 1}" for i in range(ELEMENT_COUNT)]
+    for i in range(ELEMENT_COUNT):
+        plot_single(
+            axes_masses[i],
+            axes_ratios[i],
+            f"element_mass_sp{i}",
+            columns[i],
+            **plotting_kwargs,
+        )
+
+    # Add Colorbars
+    cb_masses = fig_masses.colorbar(
+        ScalarMappable(**imshow_mass_kwargs),
+        orientation="horizontal",
+        shrink=0.75,
+        pad=0.01,
+        ax=axes_masses,
+    )
+    cb_ratios = fig_ratios.colorbar(
+        ScalarMappable(**imshow_fraction_kwargs),
+        orientation="horizontal",
+        shrink=0.75,
+        pad=0.01,
+        ax=axes_ratios,
+    )
+    cb_masses.ax.set_xlabel("Surface density (g/cm^2)")
+    cb_ratios.ax.set_xlabel("Mass ratio")
+
+    # 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/MetalAdvectionTestGEAR/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
+
+# Always generate ICs in case ELEMENT_COUNT has been changed
+echo "Generating ICs"
+python3 makeIC.py
+
+# Run SWIFT (must be compiled with --with-chemistry=GEAR_X or --with-chemistry=AGORA)
+../../../swift \
+    --hydro --threads=4 advect_metals.yml 2>&1 | tee output.log
+
+python3 runSanityChecksAdvectedMetals.py
+python3 plotAdvectedMetals.py

examples/ChemistryTests/MetalAdvectionTestGEAR/runSanityChecksAdvectedMetals.py

+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2024 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
+
+from makeIC import ELEMENT_COUNT
+
+
+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("total metal mass conservation")
+def test_total_metal_mass_conservation(data_start, data_end):
+    def metal_mass(data):
+        columns = getattr(data.gas.metal_mass_fractions, "named_columns", None)
+        if columns is not None:
+            return sum(
+                data.gas.masses * getattr(data.gas.metal_mass_fractions, c)
+                for c in columns
+            )
+        else:
+            return data.gas.masses.reshape(-1, 1) * 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_idx):
+    columns = getattr(data.gas.metal_mass_fractions, "named_columns", None)
+    if columns is not None:
+        return data.gas.masses * getattr(
+            data.gas.metal_mass_fractions, columns[element_idx]
+        )
+    else:
+        return data.gas.masses * data.gas.metal_mass_fractions[:, element_idx]
+
+
+@test("element-wise mass conservation")
+def test_element_wise_mass_conservation(data_start, data_end):
+    for i in range(ELEMENT_COUNT):
+        assert approx_equal(
+            np.sum(element_mass(data_start, i)), np.sum(element_mass(data_end, i))
+        )
+
+
+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_total_metal_mass_conservation(start, end)
+    test_element_wise_mass_conservation(start, end)
+
+    print("Done!")

examples/Cooling/ConstantCosmoTempEvolution/getGlass.sh

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

examples/Cooling/CoolingBox/README

+Runs a uniform box exposed to constant cooling rates computed with 
+the Grackle library.
+
+To run with Grackle mode 0, configure swift with:
+  ./configure --with-cooling=grackle_0 --with-grackle=$GRACKLE_ROOT
+
+To run with Grackle mode 1, configure swift with:
+  ./configure --with-cooling=grackle_1 --with-grackle=$GRACKLE_ROOT
+  
+To run with Grackle mode 2, configure swift with:
+  ./configure --with-cooling=grackle_2 --with-grackle=$GRACKLE_ROOT
+  
+To run with Grackle mode 3, configure swift with:
+  ./configure --with-cooling=grackle_3 --with-grackle=$GRACKLE_ROOT    

examples/Cooling/CoolingBox/coolingBox.yml

@@ -55,6 +55,7 @@ GrackleCooling:
   thermal_time_myr: 5                          # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
   self_shielding_method: 0                    # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
   self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
+  maximal_density_Hpcm3: -1                    # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter.
 
 GEARChemistry:
   initial_metallicity: 0.01295

examples/Cooling/CoolingBox/getGlass.sh

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

examples/Cooling/CoolingBox/run.sh

@@ -17,7 +17,7 @@ fi
 if [ ! -e CloudyData_UVB=HM2012.h5 ]
 then
     echo "Fetching the Cloudy tables required by Grackle..."
-    ../getCoolingTable.sh
+    ../getGrackleCoolingTable.sh
 fi
 
 # Run SWIFT

examples/Cooling/CoolingHeatingBox/README

+Runs a uniform box exposed to constant heating and ionization rates.
+This tests is the Iliev+2006 (ui.adsabs.harvard.edu/abs/2006MNRAS.369.1625I/) 
+Test 0 part 3, however without stopping the heating and ionization rates
+after 0.5 Myr. While the cooling is present, gas heats up and ionize.
+
+To run with GEAR-RT, configure swift with:
+  ./configure --with-cooling=grackle_1 --with-grackle=$GRACKLE_ROOT

examples/Cooling/CoolingHeatingBox/coolingBox.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841e43    # 1e10 Msol  #2.0e33     # Solar masses
+  UnitLength_in_cgs:   3.08567758e24 # Mpc        #3.0857e21  # Kiloparsecs
+  UnitVelocity_in_cgs: 1.0e5      # Kilometers 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:   5e-6  # 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:     1e-8  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            snap/snapshot # Common part of the name of output files
+  time_first:          0.         # Time of the first output (in internal units)
+  delta_time:          1e-7       # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-3 # Time between statistics output
+
+Scheduler:
+  tasks_per_cell: 64
+
+# 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: 100.       # Kelvin
+  
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./coolingBox.hdf5     # The file to read
+  periodic:   1
+  
+# Dimensionless pre-factor for the time-step condition
+LambdaCooling:
+  lambda_nH2_cgs:              1e-22 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3])
+  cooling_tstep_mult:          1.0        # Dimensionless pre-factor for the time-step condition
+
+# Cooling with Grackle 3.0
+GrackleCooling:
+  cloudy_table: CloudyData_UVB=HM2012.h5       # Name of the Cloudy Table (available on the grackle bitbucket repository)
+  with_UV_background: 0                        # Enable or not the UV background
+  redshift: 0                                  # Redshift to use (-1 means time based redshift)
+  with_metal_cooling: 0                        # Enable or not the metal cooling
+  max_steps: 10000                             # (optional) Max number of step when computing the initial composition
+  convergence_limit: 1e-2                      # (optional) Convergence threshold (relative) for initial composition
+  thermal_time_myr: 5                          # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
+  self_shielding_method: 0                     # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
+  self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
+  HydrogenFractionByMass : 1.                  # Hydrogen fraction by mass (default is 0.76)
+  use_radiative_transfer : 1                   # Arrays of ionization and heating rates are provided
+  RT_heating_rate_cgs        : 1.65117e-17     # heating         rate in units of / nHI_cgs (corresponding to a flux of 10^12 photons s−1 cm−2 , with a 10^5 K blackbody spectrum)
+  RT_HI_ionization_rate_cgs  : 1.63054e-06     # HI ionization   rate in cgs [1/s]          (corresponding to a flux of 10^12 photons s−1 cm−2 , with a 10^5 K blackbody spectrum)
+  RT_HeI_ionization_rate_cgs : 0               # HeI ionization  rate in cgs [1/s]
+  RT_HeII_ionization_rate_cgs: 0               # HeII ionization rate in cgs [1/s]
+  RT_H2_dissociation_rate_cgs: 0               # H2 dissociation rate in cgs [1/s]
+  maximal_density_Hpcm3:  -1                   # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter.
+  
+
+GEARChemistry:
+  initial_metallicity: 0.01295
+
+
+
+GEARPressureFloor:
+  jeans_factor: 0.       # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+

examples/Cooling/CoolingHeatingBox/getGlass.sh

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

examples/Cooling/CoolingHeatingBox/getResults.sh

+#!/bin/bash
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/CoolingHeatingBox_results.txt

examples/Planetary/SquareTest_2D/makeICDifferentMasses.py → examples/Cooling/CoolingHeatingBox/makeIC.py

 ###############################################################################
 # This file is part of SWIFT.
-# Copyright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
-#               2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+# Copyright (c) 2023 Yves Revaz (yves.revaz@epfl.ch)
+#                    Loic Hausammann (loic.hausammann@id.ethz.ch)
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published
@@ -19,93 +19,104 @@
 ##############################################################################
 
 import h5py
-from numpy import *
+import numpy as np
 
-# Generates a swift IC file for the Square test in a periodic box
+# Generates a SWIFT IC file with a constant density and pressure
 
 # Parameters
-L = 2 * 64  # Number of particles on the side
-gamma = 5.0 / 3.0  # Gas adiabatic index
-rho0 = 4.0  # Gas central density
-rho1 = 1.0  # Gas outskirt density
-P0 = 2.5  # Gas central pressure
-P1 = 2.5  # Gas central pressure
-vx = 0.0  # Random velocity for all particles
-vy = 0.0
-fileOutputName = "square.hdf5"
-# ---------------------------------------------------
+periodic = 1  # 1 For periodic box
+boxSize = 1  # 1 kiloparsec
+# rho = 3.2e3  # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
 
-vol = 1.0
+# Iliev test values
+rho = (
+    2471402.49842514
+)  # Density in code units (1-hydrogen atoms per cm^3 assuming Hydrogen only)
+T = 100  # Initial Temperature
 
-numPart = L * L
 
-pos_x = arange(0, 1, 1.0 / L)
-xv, yv = meshgrid(pos_x, pos_x)
-pos = zeros((numPart, 3), dtype=float)
-pos[:, 0] = xv.flatten()
-pos[:, 1] = yv.flatten()
-
-# Now we can get 2d masks!
-inside = logical_and.reduce([xv < 0.75, xv > 0.25, yv < 0.75, yv > 0.25])
+gamma = 5.0 / 3.0  # Gas adiabatic index
+fileName = "coolingBox.hdf5"
+# ---------------------------------------------------
 
-mass_in = rho0 / numPart
-mass_out = rho1 / numPart
+# defines some constants
+# need to be changed in plotResults.py too
+h_frac = 1.0  # 0.76
+mu = 4.0 / (1.0 + 3.0 * h_frac)
 
-m = ones_like(xv) * mass_out
-m[inside] = mass_in
-m = m.flatten()
+m_h_cgs = 1.67262192369e-24
+k_b_cgs = 1.380649e-16
 
-h = ones_like(m) / L
+# defines units
+unit_length = 3.08567758e24  # Mpc
+unit_mass = 1.98841e43  # 1e10 solar mass
+unit_velocity = 1e5  # 1e5 km/s
+unit_time = unit_length / unit_velocity
 
-u_in = P0 / ((gamma - 1) * rho0)
-u_out = P1 / ((gamma - 1) * rho1)
-u = ones_like(xv) * u_out
-u[inside] = u_in
-u = u.flatten()
-vel = zeros_like(pos)
-ids = arange(numPart)
-mat = zeros(numPart)
+# Read id, position and h from glass
+glass = h5py.File("glassCube_32.hdf5", "r")
+ids = glass["/PartType0/ParticleIDs"][:]
+pos = glass["/PartType0/Coordinates"][:, :] * boxSize
+h = glass["/PartType0/SmoothingLength"][:] * boxSize
 
+# Compute basic properties
+numPart = np.size(pos) // 3
+mass = boxSize ** 3 * rho / numPart
+internalEnergy = k_b_cgs * T * mu / ((gamma - 1.0) * m_h_cgs)
+internalEnergy *= (unit_time / unit_length) ** 2
 
 # File
-fileOutput = h5py.File(fileOutputName, "w")
+f = h5py.File(fileName, "w")
 
 # Header
-grp = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = [vol, vol, 0.2]
+grp = f.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
 grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [numPart, 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, 0, 0, 0, 0, 0]
-grp.attrs["Dimension"] = 2
+grp.attrs["Flag_Entropy_ICs"] = 0
+
+# Runtime parameters
+grp = f.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
 
 # Units
-grp = fileOutput.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 = f.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = unit_length
+grp.attrs["Unit mass in cgs (U_M)"] = unit_mass
+grp.attrs["Unit time in cgs (U_t)"] = unit_time
 grp.attrs["Unit current in cgs (U_I)"] = 1.0
 grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
 # Particle group
-grp = fileOutput.create_group("/PartType0")
-ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
-ds[()] = pos
+grp = f.create_group("/PartType0")
+
+v = np.zeros((numPart, 3))
 ds = grp.create_dataset("Velocities", (numPart, 3), "f")
-ds[()] = vel
+ds[()] = v
+
+m = np.full((numPart, 1), mass)
 ds = grp.create_dataset("Masses", (numPart, 1), "f")
-ds[()] = m.reshape((numPart, 1))
-ds = grp.create_dataset("SmoothingLengths", (numPart, 1), "f")
-ds[()] = h.reshape((numPart, 1))
-ds = grp.create_dataset("InternalEnergies", (numPart, 1), "f")
-ds[()] = u.reshape((numPart, 1))
+ds[()] = m
+
+h = np.reshape(h, (numPart, 1))
+ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+ds[()] = h
+
+u = np.full((numPart, 1), internalEnergy)
+ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+ds[()] = u
+
+ids = np.reshape(ids, (numPart, 1))
 ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
-ds[()] = ids.reshape((numPart, 1))
-ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
-ds[()] = mat.reshape((numPart, 1))
+ds[()] = ids
+
+ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
+ds[()] = pos
 
+f.close()
 
-fileOutput.close()
+print("Initial condition generated")

examples/Cooling/CoolingHeatingBox/plotResults.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Yves Revaz (yves.revaz@epfl.ch)
+#                    Loic Hausammann (loic.hausammann@id.ethz.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+import matplotlib
+
+import matplotlib.pyplot as plt
+import numpy as np
+from glob import glob
+import h5py
+
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
+
+
+##########################################
+# read specific energy from the solution
+##########################################
+
+# Read in units.
+resultfile = "CoolingHeatingBox_results.txt"
+
+f = open(resultfile, "r")
+firstline = f.readline()
+massline = f.readline()
+lengthline = f.readline()
+velline = f.readline()
+f.close()
+units = []
+for l in [massline, lengthline, velline]:
+    before, after = l.split("used:")
+    val, unit = after.split("[")
+    val = val.strip()
+    units.append(float(val))
+
+mass_units = units[0]
+length_units = units[1]
+velocity_units = units[2]
+time_units = velocity_units / length_units
+density_units = mass_units / length_units ** 3
+
+# Read in all other data
+data = np.loadtxt(resultfile)
+
+Time = data[:, 1]
+Time_Myr = Time * 1e-6
+Energy = data[:, 12]  # internal energy IU
+
+
+##########################################
+# compute specific energy from the models
+##########################################
+
+
+# First snapshot
+snap_filename = "snap/snapshot_0000.hdf5"
+
+# Read the initial state of the gas
+f = h5py.File(snap_filename, "r")
+
+# Read the units parameters from the snapshot
+units = f["InternalCodeUnits"]
+unit_mass = units.attrs["Unit mass in cgs (U_M)"]
+unit_length = units.attrs["Unit length in cgs (U_L)"]
+unit_time = units.attrs["Unit time in cgs (U_t)"]
+
+
+# Read snapshots
+files = glob("snap/snapshot_*.hdf5")
+files.sort()
+N = len(files)
+energy_snap = np.zeros(N)
+time_snap_cgs = np.zeros(N)
+for i in range(N):
+    snap = h5py.File(files[i], "r")
+    masses = snap["/PartType0/Masses"][:]
+    u = snap["/PartType0/InternalEnergies"][:] * masses
+    u = sum(u) / masses.sum()
+    energy_snap[i] = u
+    time_snap_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
+
+
+##########################################
+# plot
+##########################################
+
+plt.figure()
+
+Myr_in_s = 1e6 * 365.25 * 24.0 * 60.0 * 60.0
+yr_in_s = 365.25 * 24.0 * 60.0 * 60.0
+
+plt.plot(Time_Myr, Energy, "r", ms=3, label="Expected solution")
+
+plt.plot(time_snap_cgs / Myr_in_s, energy_snap, "k", ms=2)
+
+
+plt.legend(loc="right", fontsize=8, frameon=False, handlelength=3, ncol=1)
+plt.xlabel("${\\rm{Time~[Myr]}}$", labelpad=0)
+plt.ylabel(r"$\rm{Internal Energy\,\,[IU]}$")
+
+
+plt.tight_layout()
+
+plt.show()

examples/Cooling/CoolingHeatingBox/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassCube_32.hdf5 ]
+then
+    echo "Fetching initial glass file for the cooling box example..."
+    ./getGlass.sh
+fi
+if [ ! -e coolingBox.hdf5 ]
+then
+    echo "Generating initial conditions for the cooling box example..."
+    python3 makeIC.py
+fi
+
+# Get the Grackle cooling table
+if [ ! -e CloudyData_UVB=HM2012.h5 ]
+then
+    echo "Fetching the Cloudy tables required by Grackle..."
+    ../getGrackleCoolingTable.sh
+fi
+
+# Get the results
+if [ ! -e CoolingHeatingBox_results.txt]
+then
+    echo "Fetching the the results..."
+    ./getResults.sh
+fi
+
+
+# Create output directory
+rm -rf snap
+mkdir snap
+
+# Run SWIFT
+../../../swift --hydro --cooling --threads=14  coolingBox.yml
+
+# Check energy conservation and cooling rate
+python3 plotResults.py

examples/Cooling/CoolingRedshiftDependence/getGlass.sh

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

examples/Cooling/CoolingSedovBlast_3D/getGlass.sh

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