diff --git a/examples/RadiativeTransferTests/CosmoCoolingTest/README b/examples/RadiativeTransferTests/CosmoCoolingTest/README
new file mode 100644
index 0000000000000000000000000000000000000000..eb50535a856ff2a0a09f1440272ee9f22601d7e5
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoCoolingTest/README
@@ -0,0 +1,7 @@
+Check on a simple case whether the cooling works correctly without any
+radiation being present: Use a uniform 3D box of hot gas and let it cool down.
+
+The reference solution assumes case A recombination.
+
+To run with GEAR-RT, compile swift with:
+ --with-rt=GEAR_3 --with-rt-riemann-solver=GLF --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
diff --git a/examples/RadiativeTransferTests/CosmoCoolingTest/getReference.sh b/examples/RadiativeTransferTests/CosmoCoolingTest/getReference.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a6c154ac1594e94374f468f3059c6e26fb58db3d
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoCoolingTest/getReference.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/CosmoRTCoolingTestReference.txt
diff --git a/examples/RadiativeTransferTests/CosmoCoolingTest/makeIC.py b/examples/RadiativeTransferTests/CosmoCoolingTest/makeIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..b727376480e8267bd36dd15e6c0eb044a921814d
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoCoolingTest/makeIC.py
@@ -0,0 +1,158 @@
+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+# 2024 Stan Verhoeve (s06verhoeve@gmail.com)
+# 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/>.
+#
+##############################################################################
+
+# -----------------------------------------------------------
+# Use 10 particles in each dimension to generate a uniform
+# box with high temperatures.
+# -----------------------------------------------------------
+
+import h5py
+import numpy as np
+import unyt
+from swiftsimio import Writer
+from swiftsimio.units import cosmo_units
+import yaml
+
+with open(r"rt_cooling_test.yml") as paramfile:
+ params = yaml.load(paramfile, Loader=yaml.FullLoader)
+ a_begin = params["Cosmology"]["a_begin"]
+ a_begin = float(a_begin)
+
+
+# number of particles in each dimension
+n_p = 10
+nparts = n_p ** 3
+# filename of ICs to be generated
+outputfilename = "cooling_test.hdf5"
+# adiabatic index
+gamma = 5.0 / 3.0
+# total hydrogen mass fraction
+XH = 0.76
+# total helium mass fraction
+XHe = 0.24
+# boxsize
+boxsize = 1 * unyt.kpc
+# initial gas temperature
+initial_temperature = 1e6 * unyt.K
+
+# Initial particle density and mass
+gas_density_phys_cgs = 1.6756058890024518e-25 * unyt.g / unyt.cm ** 3
+
+# Include a^3 to convert physical density to comoving
+pmass = (gas_density_phys_cgs) * (boxsize ** 3 / nparts) * a_begin ** 3
+pmass = pmass.to("Msun")
+# -----------------------------------------------
+
+
+def internal_energy(T, mu):
+ """
+ Compute the internal energy of the gas for a given
+ temperature and mean molecular weight
+ """
+ # Using u = 1 / (gamma - 1) * p / rho
+ # and p = N/V * kT = rho / (mu * m_u) * kT
+
+ u = unyt.boltzmann_constant * T / (gamma - 1) / (mu * unyt.atomic_mass_unit)
+ return u
+
+
+def mean_molecular_weight(XH0, XHp, XHe0, XHep, XHepp):
+ """
+ Determines the mean molecular weight for given
+ mass fractions of
+ hydrogen: XH0
+ H+: XHp
+ He: XHe0
+ He+: XHep
+ He++: XHepp
+
+ returns:
+ mu: mean molecular weight [in atomic mass units]
+ NOTE: to get the actual mean mass, you still need
+ to multiply it by m_u, as is tradition in the formulae
+ """
+
+ # 1/mu = sum_j X_j / A_j * (1 + E_j)
+ # A_H = 1, E_H = 0
+ # A_Hp = 1, E_Hp = 1
+ # A_He = 4, E_He = 0
+ # A_Hep = 4, E_Hep = 1
+ # A_Hepp = 4, E_Hepp = 2
+ one_over_mu = XH0 + 2 * XHp + 0.25 * XHe0 + 0.5 * XHep + 0.75 * XHepp
+
+ return 1.0 / one_over_mu
+
+
+# assume everything is ionized initially
+mu = mean_molecular_weight(0.0, XH, 0.0, 0.0, XHe)
+u_part = internal_energy(initial_temperature, mu)
+pmass = pmass.to("Msun")
+
+
+xp = unyt.unyt_array(np.zeros((nparts, 3), dtype=np.float32), boxsize.units)
+dx = boxsize / n_p
+ind = 0
+for i in range(n_p):
+ x = (i + 0.5) * dx
+ for j in range(n_p):
+ y = (j + 0.5) * dx
+ for k in range(n_p):
+ z = (k + 0.5) * dx
+
+ xp[ind] = (x, y, z)
+ ind += 1
+
+w = Writer(cosmo_units, boxsize, dimension=3)
+
+w.gas.coordinates = xp
+w.gas.velocities = np.zeros(xp.shape, dtype=np.float32) * (unyt.km / unyt.s)
+w.gas.masses = np.ones(nparts, dtype=np.float32) * pmass
+w.gas.internal_energy = np.ones(nparts, dtype=np.float32) * u_part
+
+# Generate initial guess for smoothing lengths based on MIPS
+w.gas.generate_smoothing_lengths(boxsize=boxsize, dimension=3)
+
+# If IDs are not present, this automatically generates
+w.write(outputfilename)
+
+
+# Now open file back up again and add RT data.
+F = h5py.File(outputfilename, "r+")
+header = F["Header"]
+parts = F["/PartType0"]
+
+# Create initial ionization species mass fractions.
+# Assume everything is ionized initially
+# NOTE: grackle doesn't really like exact zeroes, so
+# use something very small instead.
+HIdata = np.ones(nparts, dtype=np.float32) * 1e-12
+HIIdata = np.ones(nparts, dtype=np.float32) * XH
+HeIdata = np.ones(nparts, dtype=np.float32) * 1e-12
+HeIIdata = np.ones(nparts, dtype=np.float32) * 1e-12
+HeIIIdata = np.ones(nparts, dtype=np.float32) * XHe
+
+parts.create_dataset("MassFractionHI", data=HIdata)
+parts.create_dataset("MassFractionHII", data=HIIdata)
+parts.create_dataset("MassFractionHeI", data=HeIdata)
+parts.create_dataset("MassFractionHeII", data=HeIIdata)
+parts.create_dataset("MassFractionHeIII", data=HeIIIdata)
+
+# close up, and we're done!
+F.close()
diff --git a/examples/RadiativeTransferTests/CosmoCoolingTest/plotSolution.py b/examples/RadiativeTransferTests/CosmoCoolingTest/plotSolution.py
new file mode 100755
index 0000000000000000000000000000000000000000..fb4e103ef12d1c68ebff3c6a23874ba086d565a1
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoCoolingTest/plotSolution.py
@@ -0,0 +1,475 @@
+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+# 2024 Stan Verhoeve (s06verhoeve@gmail.com)
+# 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/>.
+#
+##############################################################################
+
+# -------------------------------------------
+# Plot the gas temperature, mean molecular
+# weight, and mass fractions
+# -------------------------------------------
+
+import copy
+import os
+
+import numpy as np
+import swiftsimio
+import unyt
+from matplotlib import pyplot as plt
+
+# arguments for plots of results
+plotkwargs = {"alpha": 0.5}
+# arguments for plot of references
+referenceplotkwargs = {"color": "grey", "lw": 4, "alpha": 0.6}
+# arguments for legends
+legendprops = {"size": 8}
+# snapshot basenames
+snapshot_base = "output"
+# Reference file name
+reference_file = "CosmoRTCoolingTestReference.txt"
+# Plot time on x-axis
+plot_time = False
+
+
+# -----------------------------------------------------------------------
+energy_units = unyt.Msun * unyt.kpc ** 2 / unyt.kyr ** 2
+mass_units = unyt.Msun
+length_units = unyt.kpc
+
+
+def mean_molecular_weight(XH0, XHp, XHe0, XHep, XHepp):
+ """
+ Determines the mean molecular weight for given
+ mass fractions of
+ hydrogen: XH0
+ H+: XHp
+ He: XHe0
+ He+: XHep
+ He++: XHepp
+
+ returns:
+ mu: mean molecular weight [in atomic mass units]
+ NOTE: to get the actual mean mass, you still need
+ to multiply it by m_u, as is tradition in the formulae
+ """
+
+ # 1/mu = sum_j X_j / A_j * (1 + E_j)
+ # A_H = 1, E_H = 0
+ # A_Hp = 1, E_Hp = 1
+ # A_He = 4, E_He = 0
+ # A_Hep = 4, E_Hep = 1
+ # A_Hepp = 4, E_Hepp = 2
+ one_over_mu = XH0 + 2 * XHp + 0.25 * XHe0 + 0.5 * XHep + 0.75 * XHepp
+
+ return 1.0 / one_over_mu
+
+
+def gas_temperature(u, mu, gamma):
+ """
+ Compute the gas temperature given the specific internal
+ energy u and the mean molecular weight mu
+ """
+
+ # Using u = 1 / (gamma - 1) * p / rho
+ # and p = N/V * kT = rho / (mu * m_u) * kT
+
+ T = u * (gamma - 1) * mu * unyt.atomic_mass_unit / unyt.boltzmann_constant
+
+ return T
+
+
+def get_snapshot_list(snapshot_basename="output"):
+ """
+ Find the snapshot(s) that are to be plotted
+ and return their names as list
+ """
+
+ snaplist = []
+
+ dirlist = os.listdir()
+ for f in dirlist:
+ if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+ snaplist.append(f)
+
+ if len(snaplist) == 0:
+ raise FileNotFoundError(
+ "Didn't find any snapshots with basename '" + snapshot_basename + "'"
+ )
+
+ snaplist = sorted(snaplist)
+
+ return snaplist
+
+
+def get_ion_mass_fractions(swiftsimio_loaded_data):
+ """
+ Returns the ion mass fractions according to
+ the used scheme.
+
+ swiftsimio_loaded_data: the swiftsimio.load() object
+ """
+
+ data = swiftsimio_loaded_data
+ meta = data.metadata
+ gas = data.gas
+ with_rt = True
+ scheme = None
+ try:
+ scheme = str(meta.subgrid_scheme["RT Scheme"].decode("utf-8"))
+ except KeyError:
+ # allow to read in solutions with only cooling, without RT
+ with_rt = False
+
+ if with_rt:
+ if scheme.startswith("GEAR M1closure"):
+ imf = data.gas.ion_mass_fractions
+ elif scheme.startswith("SPH M1closure"):
+ # atomic mass
+ mamu = {
+ "e": 0.0,
+ "HI": 1.0,
+ "HII": 1.0,
+ "HeI": 4.0,
+ "HeII": 4.0,
+ "HeIII": 4.0,
+ }
+ mass_function_hydrogen = data.gas.rt_element_mass_fractions.hydrogen
+ imf = copy.deepcopy(data.gas.rt_species_abundances)
+ named_columns = data.gas.rt_species_abundances.named_columns
+ for column in named_columns:
+ # abundance is in n_X/n_H unit. We convert it to mass fraction by multipling mass fraction of H
+ mass_function = (
+ getattr(data.gas.rt_species_abundances, column)
+ * mass_function_hydrogen
+ * mamu[column]
+ )
+ setattr(imf, column, mass_function)
+ else:
+ raise ValueError("Unknown scheme", scheme)
+ else:
+ # try to find solutions for cooling only runs
+ imf = {
+ "HI": gas.hi[:],
+ "HII": gas.hii[:],
+ "HeI": gas.he_i[:],
+ "HeII": gas.he_ii[:],
+ "HeIII": gas.he_iii[:],
+ }
+
+ return imf
+
+
+def get_snapshot_data(snaplist):
+ """
+ Extract the relevant data from the list of snapshots.
+
+ Returns:
+ numpy arrays of:
+ time
+ temperatures
+ mean molecular weights
+ mass fractions
+ """
+
+ nsnaps = len(snaplist)
+ firstdata = swiftsimio.load(snaplist[0])
+ with_rt = True
+ try:
+ ngroups = int(firstdata.metadata.subgrid_scheme["PhotonGroupNumber"])
+ except KeyError:
+ # allow to read in solutions with only cooling, without RT
+ with_rt = False
+
+ # Read internal units
+ unit_system = firstdata.metadata.units
+ mass_units = unit_system.mass
+ length_units = unit_system.length
+ time_units = unit_system.time
+
+ # Units derived from base
+ velocity_units = length_units / time_units
+ energy_units = mass_units * velocity_units ** 2
+ density_units = mass_units / length_units ** 3
+
+ scale_factors = np.zeros(nsnaps)
+ times = np.zeros(nsnaps) * unyt.Myr
+ temperatures = np.zeros(nsnaps) * unyt.K
+ mean_molecular_weights = np.zeros(nsnaps)
+ mass_fractions = np.zeros((nsnaps, 5))
+ internal_energies = np.zeros(nsnaps) * energy_units
+ specific_internal_energies = np.zeros(nsnaps) * energy_units / mass_units
+ densities = np.zeros(nsnaps) * density_units
+
+ if with_rt:
+ photon_energies = np.zeros((ngroups, nsnaps)) * energy_units
+ else:
+ photon_energies = None
+
+ for i, snap in enumerate(snaplist):
+
+ data = swiftsimio.load(snap)
+ gamma = data.gas.metadata.gas_gamma[0]
+ time = data.metadata.time.copy()
+ scale_factor = data.metadata.scale_factor.copy()
+
+ gas = data.gas
+ u = gas.internal_energies[:].to(energy_units / mass_units)
+ u.convert_to_physical()
+
+ rho = gas.densities
+ rho.convert_to_physical()
+ rho.to(density_units)
+ rho.convert_to_physical()
+
+ masses = gas.masses[:].to(mass_units)
+ masses.convert_to_physical()
+
+ imf = get_ion_mass_fractions(data)
+ mu = mean_molecular_weight(imf.HI, imf.HII, imf.HeI, imf.HeII, imf.HeIII)
+ mu.convert_to_physical()
+
+ T = gas_temperature(u, mu, gamma).to("K")
+ um = u.to(energy_units / mass_units) * masses
+ um.to(energy_units)
+
+ scale_factors[i] = scale_factor
+ times[i] = time.to("Myr")
+ temperatures[i] = np.mean(T)
+ mean_molecular_weights[i] = np.mean(mu)
+ internal_energies[i] = np.mean(um)
+ specific_internal_energies[i] = u.to(energy_units / mass_units).sum() / len(u)
+ densities[i] = rho.sum() / len(rho)
+
+ mass_fractions[i, 0] = np.mean(imf.HI)
+ mass_fractions[i, 1] = np.mean(imf.HII)
+ mass_fractions[i, 2] = np.mean(imf.HeI)
+ mass_fractions[i, 3] = np.mean(imf.HeII)
+ mass_fractions[i, 4] = np.mean(imf.HeIII)
+
+ if with_rt:
+ for g in range(ngroups):
+ en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+ en = en[:].to_physical().to(energy_units)
+ photon_energies[g, i] = en.sum() / en.shape[0]
+
+ return (
+ scale_factors,
+ times,
+ temperatures,
+ mean_molecular_weights,
+ mass_fractions,
+ internal_energies,
+ photon_energies,
+ specific_internal_energies,
+ densities,
+ )
+
+
+def get_reference_data(reference_file="CosmoRTCoolingTestReference.txt"):
+ # Grab lines containing units
+ with open(reference_file) as f:
+ # Skip first line
+ f.readline()
+ mass_units_line = f.readline()
+ length_units_line = f.readline()
+ velocity_units_line = f.readline()
+
+ # Extract numbers from lines
+ mass_units = float(mass_units_line[18:-4]) * unyt.g
+ length_units = float(length_units_line[20:-5]) * unyt.cm
+ velocity_units = float(velocity_units_line[22:-7]) * unyt.cm / unyt.s
+
+ # Derived units
+ energy_units = mass_units * velocity_units ** 2
+ density_units = mass_units / length_units ** 3
+
+ # Read in data
+ refdata = np.loadtxt(reference_file, dtype=np.float64)
+
+ a_ref = refdata[:, 1]
+ t_ref = refdata[:, 3] * unyt.yr
+ dt_ref = refdata[:, 4] * unyt.yr
+ T_ref = refdata[:, 5] * unyt.K
+ mu_ref = refdata[:, 6]
+ rho_ref = refdata[:, 7] * density_units
+ rhoHI_ref = refdata[:, 8] * density_units
+ rhoHII_ref = refdata[:, 9] * density_units
+ rhoHeI_ref = refdata[:, 10] * density_units
+ rhoHeII_ref = refdata[:, 11] * density_units
+ rhoHeIII_ref = refdata[:, 12] * density_units
+ rhoe_ref = refdata[:, 13] * density_units
+ u_ref = refdata[:, -1] * energy_units / mass_units
+ mass_fraction_ref = np.empty((t_ref.shape[0], 5))
+ mass_fraction_ref[:, 0] = rhoHI_ref / rho_ref
+ mass_fraction_ref[:, 1] = rhoHII_ref / rho_ref
+ mass_fraction_ref[:, 2] = rhoHeI_ref / rho_ref
+ mass_fraction_ref[:, 3] = rhoHeII_ref / rho_ref
+ mass_fraction_ref[:, 4] = rhoHeIII_ref / rho_ref
+
+ return (
+ a_ref,
+ t_ref,
+ T_ref,
+ mu_ref,
+ mass_fraction_ref,
+ dt_ref,
+ rho_ref,
+ rhoHI_ref,
+ rhoHII_ref,
+ rhoHeI_ref,
+ rhoHeII_ref,
+ rhoHeIII_ref,
+ rhoe_ref,
+ u_ref,
+ )
+
+
+if __name__ == "__main__":
+
+ # Get list of shapshots
+ snaplist = get_snapshot_list(snapshot_base)
+
+ # Read snapshot data
+ a, t, T, mu, mass_fraction, u, photon_energies, us, rho = get_snapshot_data(
+ snaplist
+ )
+
+ with_rt = photon_energies is not None
+ if with_rt:
+ ngroups = photon_energies.shape[0]
+
+ # Read reference solution data
+ a_ref, t_ref, T_ref, mu_ref, mass_fraction_ref, *__ = get_reference_data(
+ reference_file
+ )
+
+ # Convert t_ref to Myr
+ t_ref.convert_to_units("Myr")
+
+ # Translate snapshot times to start at t=0
+ t -= t[0]
+
+ # Grab x-coordinate
+ if plot_time:
+ xcoords = t
+ xcoords_ref = t_ref
+ xlabel = "Time [Myr]"
+ xscale = "log"
+ xlims = (0.1, max(t))
+ else:
+ xcoords = a
+ xcoords_ref = a_ref
+ xlabel = "Scale factor [1]"
+ xscale = "linear"
+ xlims = (min(a), max(a))
+
+ # ------------------
+ # Plot figures
+ # ------------------
+
+ fig = plt.figure(figsize=(8, 8), dpi=300)
+ ax1 = fig.add_subplot(2, 2, 1)
+ ax2 = fig.add_subplot(2, 2, 2)
+ ax3 = fig.add_subplot(2, 2, 3)
+ ax4 = fig.add_subplot(2, 2, 4)
+
+ ax1.set_xscale(xscale)
+ ax2.set_xscale(xscale)
+ ax3.set_xscale(xscale)
+ ax4.set_xscale(xscale)
+
+ ax1.set_xlim(*xlims)
+ ax2.set_xlim(*xlims)
+ ax3.set_xlim(*xlims)
+ ax4.set_xlim(*xlims)
+
+ ax1.semilogy(xcoords_ref, T_ref, label="reference", **referenceplotkwargs)
+ ax1.semilogy(xcoords, T, label="obtained results")
+
+ ax1.set_yscale("log")
+ ax1.set_xlabel(xlabel)
+ ax1.set_ylabel("gas temperature [K]")
+ ax1.legend(prop=legendprops)
+ ax1.grid()
+
+ ax2.plot(xcoords_ref, mu_ref, label="reference", **referenceplotkwargs)
+ ax2.plot(xcoords, mu, label="obtained results")
+
+ ax2.set_xlabel(xlabel)
+ ax2.set_ylabel("mean molecular weight")
+ ax2.legend(prop=legendprops)
+ ax2.grid()
+
+ total_mass_fraction = np.sum(mass_fraction, axis=1)
+ ax3.plot(xcoords, total_mass_fraction, "k", label="total", ls="-")
+
+ ax3.plot(
+ xcoords_ref[1:],
+ mass_fraction_ref[1:, 0],
+ label="reference",
+ **referenceplotkwargs,
+ zorder=0,
+ )
+ ax3.plot(xcoords_ref[1:], mass_fraction_ref[1:, 1], **referenceplotkwargs, zorder=0)
+ ax3.plot(xcoords_ref[1:], mass_fraction_ref[1:, 2], **referenceplotkwargs, zorder=0)
+ ax3.plot(xcoords_ref[1:], mass_fraction_ref[1:, 3], **referenceplotkwargs, zorder=0)
+ ax3.plot(xcoords_ref[1:], mass_fraction_ref[1:, 4], **referenceplotkwargs, zorder=0)
+
+ ax3.plot(xcoords, mass_fraction[:, 0], label="HI", ls=":", **plotkwargs, zorder=1)
+ ax3.plot(xcoords, mass_fraction[:, 1], label="HII", ls="-.", **plotkwargs, zorder=1)
+ ax3.plot(xcoords, mass_fraction[:, 2], label="HeI", ls=":", **plotkwargs, zorder=1)
+ ax3.plot(
+ xcoords, mass_fraction[:, 3], label="HeII", ls="-.", **plotkwargs, zorder=1
+ )
+ ax3.plot(
+ xcoords, mass_fraction[:, 4], label="HeIII", ls="--", **plotkwargs, zorder=1
+ )
+ ax3.legend(loc="upper right", prop=legendprops)
+ ax3.set_xlabel(xlabel)
+ ax3.set_ylabel("gas mass fractions [1]")
+ ax3.grid()
+
+ if with_rt:
+ u.convert_to_units(energy_units)
+ photon_energies.convert_to_units(energy_units)
+ tot_energy = u + np.sum(photon_energies, axis=0)
+ ax4.plot(
+ xcoords,
+ tot_energy,
+ label=f"total energy budget",
+ color="k",
+ ls="--",
+ **plotkwargs,
+ )
+ for g in range(ngroups):
+ ax4.plot(
+ xcoords,
+ photon_energies[g, :],
+ label=f"photon energies group {g + 1}",
+ **plotkwargs,
+ )
+ ax4.plot(xcoords, u, label="gas internal energy", **plotkwargs)
+ ax4.set_xlabel(xlabel)
+ ax4.set_ylabel(
+ r"energy budget [$" + u.units.latex_representation() + "$]", usetex=True
+ )
+ ax4.legend(prop=legendprops)
+ ax4.grid()
+
+ plt.tight_layout()
+ # plt.show()
+ plt.savefig("cooling_test.png")
diff --git a/examples/RadiativeTransferTests/CosmoCoolingTest/rt_cooling_test.yml b/examples/RadiativeTransferTests/CosmoCoolingTest/rt_cooling_test.yml
new file mode 100644
index 0000000000000000000000000000000000000000..97c3724cb5c51a21d36efee861ce1367ceed400f
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoCoolingTest/rt_cooling_test.yml
@@ -0,0 +1,80 @@
+MetaData:
+ run_name: RT Cooling Test
+
+# 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 # 1 kpc in cm
+ UnitVelocity_in_cgs: 1e5 # 1 km/s in cm/s
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+
+# Parameters governing the time integration
+TimeIntegration:
+ max_nr_rt_subcycles: 1
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 0.100 # The end time of the simulation (in internal units).
+ dt_min: 1.e-8 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 4.882814e-03 # 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
+ output_list_on: 0
+ output_list: snaplist_cooling
+ scale_factor_first: 0.047601 # Time of the first output (in internal units)
+ delta_time: 1.006
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ scale_factor_first: 0.047601
+ delta_time: 1.006 # 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.6 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 0. # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./cooling_test.hdf5 # The file to read
+ periodic: 1 # periodic ICs
+
+GEARRT:
+ f_reduce_c: 1.e-9 # This test is without actual radiation, so we don't care about this
+ CFL_condition: 0.9 # CFL condition for RT, independent of hydro
+ photon_groups_Hz: [3.288e15, 5.945e15, 13.157e15] # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+ stellar_luminosity_model: const # Which model to use to determine the stellar luminosities.
+ const_stellar_luminosities_LSol: [0., 0., 0.] # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+ hydrogen_mass_fraction: 0.76 # total hydrogen mass fraction
+ stellar_spectrum_type: 0 # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+ stellar_spectrum_const_max_frequency_Hz: 1.e17 # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum.
+ case_B_recombination: 0 # reference solution assumes case A recombination
+
+
+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
+ provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates
+ provide_specific_heating_rates: 0 # (optional) User provide specific heating rates
+ max_steps: 10000 # (optional) Max number of step when computing the initial composition
+ convergence_limit: 1 # (optional) Convergence threshold (relative) for initial composition
+ self_shielding_method: 0 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
+ primordial_chemistry: 1
+ thermal_time_myr: 5
+
+Scheduler:
+ tasks_per_cell: 128
+
+Cosmology: # Planck13 (EAGLE flavour)
+ a_begin: 0.0476 # z~20
+ a_end: 0.2 # z~4
+ h: 0.6777
+ Omega_cdm: 0.2587481
+ Omega_lambda: 0.693
+ Omega_b: 0.0482519
diff --git a/examples/RadiativeTransferTests/CosmoCoolingTest/run.sh b/examples/RadiativeTransferTests/CosmoCoolingTest/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..890a4b745e973c212ac5b7ddb1b13d932607c2f4
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoCoolingTest/run.sh
@@ -0,0 +1,35 @@
+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -f ./cooling_test.hdf5 ]; then
+ echo "creating ICs"
+ python3 makeIC.py
+fi
+
+# Run SWIFT with RT
+../../../swift \
+ --hydro \
+ --cosmology \
+ --threads=4 \
+ --verbose=0 \
+ --radiation \
+ --external-gravity \
+ --stars \
+ --feedback \
+./rt_cooling_test.yml 2>&1 | tee output.log
+
+# Wanna run with cooling, but no RT? This should do the trick
+# ../../../swift \
+# --hydro \
+# --threads=4 \
+# --verbose=0 \
+# --cooling \
+# ./rt_cooling_test.yml 2>&1 | tee output.log
+
+if [ ! -f "CosmoRTCoolingTestReference.txt" ]; then
+ ./getReference.sh
+fi
+python3 plotSolution.py
diff --git a/examples/RadiativeTransferTests/CosmoHeatingTest/README b/examples/RadiativeTransferTests/CosmoHeatingTest/README
new file mode 100644
index 0000000000000000000000000000000000000000..18ba3a98b493955657afb9b6cdc309af72fbb4cd
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoHeatingTest/README
@@ -0,0 +1,6 @@
+Runs a uniform box with radiation initially present. No further radiation
+sources are present, and the gas should heat up and ionize, while the
+radiation fields should decrease at the same rate.
+
+To run with GEAR-RT, compile swift with:
+ --with-rt=GEAR_3 --with-rt-riemann-solver=GLF --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
diff --git a/examples/RadiativeTransferTests/CosmoHeatingTest/makeIC.py b/examples/RadiativeTransferTests/CosmoHeatingTest/makeIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..813c0d8a4b9f66e3274a1f574212e0adbd3bd068
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoHeatingTest/makeIC.py
@@ -0,0 +1,194 @@
+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+# 2024 Stan Verhoeve (s06verhoeve@gmail.com)
+# 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/>.
+#
+##############################################################################
+
+
+# -----------------------------------------------------------
+# Use 10 particles in each dimension to generate a uniform box
+# with a fixed amount of radiation.
+# -----------------------------------------------------------
+
+from swiftsimio import Writer
+from swiftsimio.units import cosmo_units
+
+import unyt
+import numpy as np
+import h5py
+import yaml
+
+# Grab scale factor at start of run
+with open(r"rt_heating_test.yml") as paramfile:
+ params = yaml.load(paramfile, Loader=yaml.FullLoader)
+ a_begin = params["Cosmology"]["a_begin"]
+ a_begin = float(a_begin)
+
+# Reduced speed of light
+f_reduce_speed_of_light = 1.0e-6
+
+# number of particles in each dimension
+n_p = 10
+nparts = n_p ** 3
+
+# filename of ICs to be generated
+outputfilename = "heating_test.hdf5"
+# adiabatic index
+gamma = 5.0 / 3.0
+# total hydrogen mass fraction
+XH = 0.76
+# total helium mass fraction
+XHe = 0.24
+# boxsize
+boxsize = 1 * unyt.kpc
+# initial gas temperature
+initial_temperature = 1e3 * unyt.K
+# particle mass
+pmass = (unyt.atomic_mass_unit / unyt.cm ** 3) * (boxsize ** 3 / nparts) * a_begin ** 3
+pmass = pmass.to("Msun")
+
+# -----------------------------------------------
+
+
+def internal_energy(T, mu):
+ """
+ Compute the internal energy of the gas for a given
+ temperature and mean molecular weight
+ """
+ # Using u = 1 / (gamma - 1) * p / rho
+ # and p = N/V * kT = rho / (mu * m_u) * kT
+
+ u = unyt.boltzmann_constant * T / (gamma - 1) / (mu * unyt.atomic_mass_unit)
+ return u
+
+
+def mean_molecular_weight(XH0, XHp, XHe0, XHep, XHepp):
+ """
+ Determines the mean molecular weight for given
+ mass fractions of
+ hydrogen: XH0
+ H+: XHp
+ He: XHe0
+ He+: XHep
+ He++: XHepp
+
+ returns:
+ mu: mean molecular weight [in atomic mass units]
+ NOTE: to get the actual mean mass, you still need
+ to multiply it by m_u, as is tradition in the formulae
+ """
+
+ # 1/mu = sum_j X_j / A_j * (1 + E_j)
+ # A_H = 1, E_H = 0
+ # A_Hp = 1, E_Hp = 1
+ # A_He = 4, E_He = 0
+ # A_Hep = 4, E_Hep = 1
+ # A_Hepp = 4, E_Hepp = 2
+ one_over_mu = XH0 + 2 * XHp + 0.25 * XHe0 + 0.5 * XHep + 0.75 * XHepp
+
+ return 1.0 / one_over_mu
+
+
+# assume everything is neutral initially
+mu = mean_molecular_weight(XH, 0, XHe, 0.0, 0.0)
+u_part = internal_energy(initial_temperature, mu)
+pmass = pmass.to("Msun")
+
+
+xp = unyt.unyt_array(np.zeros((nparts, 3), dtype=np.float32), boxsize.units)
+dx = boxsize / n_p
+ind = 0
+for i in range(n_p):
+ x = (i + 0.5) * dx
+ for j in range(n_p):
+ y = (j + 0.5) * dx
+ for k in range(n_p):
+ z = (k + 0.5) * dx
+
+ xp[ind] = (x, y, z)
+ ind += 1
+
+w = Writer(cosmo_units, boxsize, dimension=3)
+
+
+w.gas.coordinates = xp
+w.gas.velocities = np.zeros(xp.shape, dtype=np.float32) * (unyt.km / unyt.s)
+w.gas.masses = np.ones(nparts, dtype=np.float32) * pmass
+w.gas.internal_energy = np.ones(nparts, dtype=np.float32) * u_part
+
+# Generate initial guess for smoothing lengths based on MIPS
+w.gas.generate_smoothing_lengths(boxsize=boxsize, dimension=3)
+
+# If IDs are not present, this automatically generates
+w.write(outputfilename)
+
+# Now open file back up again and add RT data.
+
+F = h5py.File(outputfilename, "r+")
+header = F["Header"]
+nparts = header.attrs["NumPart_ThisFile"][0]
+parts = F["/PartType0"]
+
+# Create initial ionization species mass fractions.
+# Assume everything neutral initially
+# NOTE: grackle doesn't really like exact zeroes, so
+# use something very small instead.
+HIdata = np.ones(nparts, dtype=np.float32) * XH
+HIIdata = np.ones(nparts, dtype=np.float32) * 1e-12
+HeIdata = np.ones(nparts, dtype=np.float32) * XHe
+HeIIdata = np.ones(nparts, dtype=np.float32) * 1e-12
+HeIIIdata = np.ones(nparts, dtype=np.float32) * 1e-12
+
+parts.create_dataset("MassFractionHI", data=HIdata)
+parts.create_dataset("MassFractionHII", data=HIIdata)
+parts.create_dataset("MassFractionHeI", data=HeIdata)
+parts.create_dataset("MassFractionHeII", data=HeIIdata)
+parts.create_dataset("MassFractionHeIII", data=HeIIIdata)
+
+
+# Add photon groups
+nPhotonGroups = 3
+
+# with this IC, the radiative cooling is negligible.
+# photon_energy = u_part * pmass * 5.0
+# photon_energy = np.arange(1, nPhotonGroups+1) * photon_energy
+
+# Fluxes from the Iliev Test0 part3
+fluxes_iliev = np.array([1.350e1, 2.779e1, 6.152e0]) * unyt.erg / unyt.s / unyt.cm ** 2
+energy_density = fluxes_iliev / unyt.c
+photon_energy = energy_density * boxsize ** 3 / nparts * a_begin ** 3
+photon_energy = photon_energy * f_reduce_speed_of_light
+
+photon_energy.convert_to_units(cosmo_units["energy"])
+photon_fluxes = 0.333333 * unyt.c * photon_energy
+photon_fluxes.convert_to_units(
+ cosmo_units["energy"] * cosmo_units["length"] / cosmo_units["time"]
+)
+
+
+for grp in range(nPhotonGroups):
+ dsetname = "PhotonEnergiesGroup{0:d}".format(grp + 1)
+ energydata = np.ones(nparts, dtype=np.float32) * photon_energy[grp]
+ parts.create_dataset(dsetname, data=energydata)
+
+ dsetname = "PhotonFluxesGroup{0:d}".format(grp + 1)
+ fluxdata = np.zeros((nparts, 3), dtype=np.float32)
+ fluxdata[:, 0] = photon_fluxes[grp]
+ parts.create_dataset(dsetname, data=fluxdata)
+
+# close up, and we're done!
+F.close()
diff --git a/examples/RadiativeTransferTests/CosmoHeatingTest/plotSolution.py b/examples/RadiativeTransferTests/CosmoHeatingTest/plotSolution.py
new file mode 100755
index 0000000000000000000000000000000000000000..c596c24946d0e8e25fb23b18024c9f5290a95155
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoHeatingTest/plotSolution.py
@@ -0,0 +1,310 @@
+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+# 2024 Stan Verhoeve (s06verhoeve@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# -------------------------------------------
+# Plot the gas temperature, mean molecular
+# weight, and mass fractions
+# -------------------------------------------
+
+import copy
+import os
+
+import numpy as np
+import swiftsimio
+import unyt
+from matplotlib import pyplot as plt
+from matplotlib import ticker as mticker
+
+# arguments for plots of results
+plotkwargs = {"alpha": 0.5}
+# arguments for plot of references
+referenceplotkwargs = {"color": "grey", "lw": 4, "alpha": 0.6}
+# arguments for legends
+legendprops = {"size": 8}
+# snapshot basenames
+snapshot_base = "output"
+# Plot time on x-axis
+plot_time = False
+
+# -----------------------------------------------------------------------
+
+energy_units = unyt.Msun * unyt.kpc ** 2 / unyt.kyr ** 2
+mass_units = unyt.Msun
+
+
+def mean_molecular_weight(XH0, XHp, XHe0, XHep, XHepp):
+ """
+ Determines the mean molecular weight for given
+ mass fractions of
+ hydrogen: XH0
+ H+: XHp
+ He: XHe0
+ He+: XHep
+ He++: XHepp
+
+ returns:
+ mu: mean molecular weight [in atomic mass units]
+ NOTE: to get the actual mean mass, you still need
+ to multiply it by m_u, as is tradition in the formulae
+ """
+
+ # 1/mu = sum_j X_j / A_j * (1 + E_j)
+ # A_H = 1, E_H = 0
+ # A_Hp = 1, E_Hp = 1
+ # A_He = 4, E_He = 0
+ # A_Hep = 4, E_Hep = 1
+ # A_Hepp = 4, E_Hepp = 2
+ one_over_mu = XH0 + 2 * XHp + 0.25 * XHe0 + 0.5 * XHep + 0.75 * XHepp
+
+ return 1.0 / one_over_mu
+
+
+def gas_temperature(u, mu, gamma):
+ """
+ Compute the gas temperature given the specific internal
+ energy u and the mean molecular weight mu
+ """
+
+ # Using u = 1 / (gamma - 1) * p / rho
+ # and p = N/V * kT = rho / (mu * m_u) * kT
+
+ T = u * (gamma - 1) * mu * unyt.atomic_mass_unit / unyt.boltzmann_constant
+
+ return T
+
+
+def get_snapshot_list(snapshot_basename="output"):
+ """
+ Find the snapshot(s) that are to be plotted
+ and return their names as list
+ """
+
+ snaplist = []
+
+ dirlist = os.listdir()
+ for f in dirlist:
+ if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+ snaplist.append(f)
+
+ if len(snaplist) == 0:
+ raise FileNotFoundError(
+ "Didn't find any snapshots with basename '" + snapshot_basename + "'"
+ )
+
+ snaplist = sorted(snaplist)
+
+ return snaplist
+
+
+def get_ion_mass_fractions(swiftsimio_loaded_data):
+ """
+ Returns the ion mass fractions according to
+ the used scheme.
+
+ swiftsimio_loaded_data: the swiftsimio.load() object
+ """
+
+ data = swiftsimio_loaded_data
+ meta = data.metadata
+ try:
+ scheme = str(meta.subgrid_scheme["RT Scheme"].decode("utf-8"))
+ except KeyError:
+ raise ValueError("This test needs to be run with RT on.")
+
+ if scheme.startswith("GEAR M1closure"):
+ imf = data.gas.ion_mass_fractions
+ elif scheme.startswith("SPH M1closure"):
+ # atomic mass
+ mamu = {"e": 0.0, "HI": 1.0, "HII": 1.0, "HeI": 4.0, "HeII": 4.0, "HeIII": 4.0}
+ mass_function_hydrogen = data.gas.rt_element_mass_fractions.hydrogen
+ imf = copy.deepcopy(data.gas.rt_species_abundances)
+ named_columns = data.gas.rt_species_abundances.named_columns
+ for column in named_columns:
+ # abundance is in n_X/n_H unit. We convert it to mass fraction by multipling mass fraction of H
+ mass_function = (
+ getattr(data.gas.rt_species_abundances, column)
+ * mass_function_hydrogen
+ * mamu[column]
+ )
+ setattr(imf, column, mass_function)
+ else:
+ raise ValueError("Unknown scheme", scheme)
+
+ return imf
+
+
+def get_snapshot_data(snaplist):
+ """
+ Extract the relevant data from the list of snapshots.
+
+ Returns:
+ numpy arrays of:
+ time
+ temperatures
+ mean molecular weights
+ mass fractions
+ """
+
+ nsnaps = len(snaplist)
+ firstdata = swiftsimio.load(snaplist[0])
+ ngroups = int(firstdata.metadata.subgrid_scheme["PhotonGroupNumber"])
+
+ scale_factors = np.zeros(nsnaps)
+ times = np.zeros(nsnaps) * unyt.Myr
+ temperatures = np.zeros(nsnaps) * unyt.K
+ mean_molecular_weights = np.zeros(nsnaps) * unyt.atomic_mass_unit
+ mass_fractions = np.zeros((nsnaps, 5))
+ internal_energies = np.zeros(nsnaps) * energy_units
+ photon_energies = np.zeros((ngroups, nsnaps)) * energy_units
+
+ time_first = firstdata.metadata.time.copy()
+ for i, snap in enumerate(snaplist):
+
+ data = swiftsimio.load(snap)
+ gamma = data.gas.metadata.gas_gamma[0]
+ time = data.metadata.time.copy()
+ scale_factor = data.metadata.scale_factor
+ gas = data.gas
+
+ u = gas.internal_energies.to(energy_units / mass_units)
+ u.convert_to_physical()
+ u_phys_cgs = u.to(unyt.erg / unyt.g)
+ masses = gas.masses.to(mass_units)
+ imf = get_ion_mass_fractions(data)
+ mu = mean_molecular_weight(imf.HI, imf.HII, imf.HeI, imf.HeII, imf.HeIII)
+ mu.convert_to_physical()
+ T = gas_temperature(u_phys_cgs, mu, gamma).to("K")
+ um = u.to(energy_units / mass_units) * masses
+
+ scale_factors[i] = scale_factor
+ times[i] = time.to("Myr")
+ temperatures[i] = np.mean(T)
+ mean_molecular_weights[i] = np.mean(mu)
+ internal_energies[i] = np.mean(um)
+
+ mass_fractions[i, 0] = np.mean(imf.HI)
+ mass_fractions[i, 1] = np.mean(imf.HII)
+ mass_fractions[i, 2] = np.mean(imf.HeI)
+ mass_fractions[i, 3] = np.mean(imf.HeII)
+ mass_fractions[i, 4] = np.mean(imf.HeIII)
+
+ for g in range(ngroups):
+ en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+ en = en.to(energy_units)
+ photon_energies[g, i] = en.sum() / en.shape[0]
+
+ return (
+ scale_factors,
+ times,
+ temperatures,
+ mean_molecular_weights,
+ mass_fractions,
+ internal_energies,
+ photon_energies,
+ )
+
+
+if __name__ == "__main__":
+ # ------------------
+ # Plot figures
+ # ------------------
+
+ snaplist = get_snapshot_list(snapshot_base)
+ a, t, T, mu, mass_fraction, u, photon_energies = get_snapshot_data(snaplist)
+ ngroups = photon_energies.shape[0]
+
+ if plot_time:
+ xcoords = t
+ xlabel = "Time [Myr]"
+ xscale = "log"
+ else:
+ xcoords = a
+ xlabel = "Scale factor [1]"
+ xscale = "linear"
+
+ fig = plt.figure(figsize=(8, 8), dpi=300)
+ ax1 = fig.add_subplot(2, 2, 1)
+ ax2 = fig.add_subplot(2, 2, 2)
+ ax3 = fig.add_subplot(2, 2, 3)
+ ax4 = fig.add_subplot(2, 2, 4)
+
+ ax1.semilogy(xcoords, T, label="obtained results")
+ ax1.set_xlabel(xlabel)
+ ax1.set_ylabel("gas temperature [K]")
+ ax1.legend(prop=legendprops)
+ ax1.grid()
+
+ ax2.plot(xcoords, mu, label="obtained results")
+ ax2.set_xlabel(xlabel)
+ ax2.set_ylabel("mean molecular weight")
+ ax2.legend(prop=legendprops)
+ ax2.grid()
+
+ total_mass_fraction = np.sum(mass_fraction, axis=1)
+ ax3.plot(xcoords, total_mass_fraction, "k", label="total", ls="-")
+
+ ax3.plot(xcoords, mass_fraction[:, 0], label="HI", ls=":", **plotkwargs, zorder=1)
+ ax3.plot(xcoords, mass_fraction[:, 1], label="HII", ls="-.", **plotkwargs, zorder=1)
+ ax3.plot(xcoords, mass_fraction[:, 2], label="HeI", ls=":", **plotkwargs, zorder=1)
+ ax3.plot(
+ xcoords, mass_fraction[:, 3], label="HeII", ls="-.", **plotkwargs, zorder=1
+ )
+ ax3.plot(
+ xcoords, mass_fraction[:, 4], label="HeIII", ls="--", **plotkwargs, zorder=1
+ )
+ ax3.legend(loc="upper right", prop=legendprops)
+ ax3.set_xlabel(xlabel)
+ ax3.set_ylabel("gas mass fractions [1]")
+ ax3.grid()
+
+ tot_energy = u + np.sum(photon_energies, axis=0)
+ ax4.plot(
+ xcoords,
+ tot_energy,
+ label=f"total energy budget",
+ color="k",
+ ls="--",
+ **plotkwargs,
+ )
+ for g in range(ngroups):
+ ax4.plot(
+ xcoords,
+ photon_energies[g, :],
+ label=f"photon energies group {g+1}",
+ **plotkwargs,
+ )
+ ax4.plot(xcoords, u, label=r"gas internal energy", **plotkwargs)
+
+ ax4.set_yscale("log")
+ ax4.set_xlabel(xlabel)
+ ax4.set_ylabel(
+ r"energy budget [$" + energy_units.units.latex_representation() + "$]",
+ usetex=True,
+ )
+ ax4.legend(prop=legendprops)
+ ax4.grid()
+
+ for ax in fig.axes:
+ ax.set_xscale(xscale)
+ ax.xaxis.set_minor_formatter(mticker.ScalarFormatter())
+
+ plt.tight_layout()
+ plt.savefig("heating_test.png")
diff --git a/examples/RadiativeTransferTests/CosmoHeatingTest/rt_heating_test.yml b/examples/RadiativeTransferTests/CosmoHeatingTest/rt_heating_test.yml
new file mode 100644
index 0000000000000000000000000000000000000000..e0b91d8b80f7e86986fbcc6d84ab8895541f8b08
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoHeatingTest/rt_heating_test.yml
@@ -0,0 +1,65 @@
+MetaData:
+ run_name: Heating Test
+
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.98848e33 # 1 M_sun in grams
+ UnitLength_in_cgs: 3.08567758e18 # 1 pc in cm
+ UnitVelocity_in_cgs: 1e5 # 1 km/s in cm/s
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+
+# Parameters governing the time integration
+TimeIntegration:
+ max_nr_rt_subcycles: 1
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ dt_min: 1.0208453377e-08 # 0.01 yr
+ dt_max: 0.0010208453 # 1 kyr
+ time_end: 4.10084 # 4 Myr
+
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: output # Common part of the name of output files
+ scale_factor_first: 0.00990099 # Time of the first output (in internal units)
+ output_list_on: 0 # (Optional) Enable the output list
+ output_list: snaplist.txt # (Optional) File containing the output times (see documentation in "Parameter File" section)
+ delta_time: 1.001
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ time_first: 0.00990099
+ delta_time: 1.001
+
+# 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.6 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 10. # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./heating_test.hdf5 # The file to read
+ periodic: 1 # periodic ICs
+
+GEARRT:
+ f_reduce_c: 1e-6 # We don't care about the radiation propagation in this test, so let's speed things up
+ CFL_condition: 0.9 # CFL condition for RT, independent of hydro
+ photon_groups_Hz: [3.288e15, 5.945e15, 13.157e15] # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+ stellar_luminosity_model: const # Which model to use to determine the stellar luminosities.
+ const_stellar_luminosities_LSol: [1., 1., 1.] # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+ set_equilibrium_initial_ionization_mass_fractions: 0 # (Optional) set the initial ionization fractions depending on gas temperature assuming ionization equilibrium.
+ hydrogen_mass_fraction: 0.76 # total hydrogen mass fraction
+ stellar_spectrum_type: 1 # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+ stellar_spectrum_blackbody_temperature_K: 1.e5 # (Conditional) if stellar_spectrum_type=1, use this temperature (in K) for the blackbody spectrum.
+ case_B_recombination: 0 # (Optional) use case B recombination interaction rates.
+
+Cosmology: # Planck13 (EAGLE flavour)
+ a_begin: 0.00990099 # z~100
+ a_end: 0.011 # z~90
+ h: 0.6777
+ Omega_cdm: 0.2587481
+ Omega_lambda: 0.693
+ Omega_b: 0.0482519
+
diff --git a/examples/RadiativeTransferTests/CosmoHeatingTest/run.sh b/examples/RadiativeTransferTests/CosmoHeatingTest/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..7a0a1e545963188c2bd3bfeaa04b3177e1de39a5
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoHeatingTest/run.sh
@@ -0,0 +1,24 @@
+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -f ./heating_test.hdf5 ]; then
+ echo "creating ICs"
+ python3 makeIC.py
+fi
+
+# Run SWIFT with RT
+../../../swift \
+ --hydro \
+ --cosmology \
+ --threads=4 \
+ --verbose=0 \
+ --radiation \
+ --external-gravity \
+ --stars \
+ --feedback \
+ ./rt_heating_test.yml 2>&1 | tee output.log
+
+python3 plotSolution.py
diff --git a/src/engine.c b/src/engine.c
index 35d60b0281a8c5b2ab5b26f670bce614a0f5bf81..f3802439661fb4200cbbea050a3303789f249e31 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1931,7 +1931,10 @@ void engine_run_rt_sub_cycles(struct engine *e) {
e->max_active_bin_subcycle = get_max_active_bin(e->ti_current_subcycle);
e->min_active_bin_subcycle =
get_min_active_bin(e->ti_current_subcycle, ti_subcycle_old);
- /* TODO: add rt_props_update() for cosmological thermochemistry*/
+
+ /* Update rt properties */
+ rt_props_update(e->rt_props, e->internal_units, e->cosmology);
+
if (e->policy & engine_policy_cosmology) {
double time_old = time;
cosmology_update(
@@ -2080,6 +2083,9 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
cooling_update(e->physical_constants, e->cosmology, e->pressure_floor_props,
e->cooling_func, e->s, e->time);
+ if (e->policy & engine_policy_rt)
+ rt_props_update(e->rt_props, e->internal_units, e->cosmology);
+
#ifdef WITH_CSDS
if (e->policy & engine_policy_csds) {
/* Mark the first time step in the particle csds file. */
@@ -2513,6 +2519,10 @@ int engine_step(struct engine *e) {
hydro_props_update(e->hydro_properties, e->gravity_properties,
e->cosmology);
+ /* Update the rt properties */
+ if (e->policy & engine_policy_rt)
+ rt_props_update(e->rt_props, e->internal_units, e->cosmology);
+
/* Check for any snapshot triggers */
engine_io_check_snapshot_triggers(e);
@@ -3995,7 +4005,7 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
struct rt_props *rt_properties =
(struct rt_props *)malloc(sizeof(struct rt_props));
rt_struct_restore(rt_properties, stream, e->physical_constants,
- e->internal_units);
+ e->internal_units, cosmo);
e->rt_props = rt_properties;
struct black_holes_props *black_holes_properties =
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 3d2e3cbded1866cf6c94c53aa86d8e8b46a76c2a..25a3d9881dd3d077e67d3d56f5f381a3bd983a8b 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -216,7 +216,7 @@ INLINE static void hydro_write_particles(const struct part* parts,
list[4] = io_make_output_field_convert_part(
"InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS,
- 3.f * hydro_gamma_minus_one, parts, xparts, convert_u,
+ -3.f * hydro_gamma_minus_one, parts, xparts, convert_u,
"Co-moving thermal energies per unit mass of the particles");
list[5] =
diff --git a/src/rt/GEAR/rt_getters.h b/src/rt/GEAR/rt_getters.h
index dc9e3e6b9cb4a5337e357aad05cb4da1976cb4bb..a11ee9cbbecc1d660160d9e2594caf8d9eea8b80 100644
--- a/src/rt/GEAR/rt_getters.h
+++ b/src/rt/GEAR/rt_getters.h
@@ -29,20 +29,35 @@
*/
/**
- * @brief Get the radiation energy densities of a particle.
+ * @brief Get the comoving radiation energy densities of a particle.
*
* @param p Particle.
* @param E (return) Pointer to the array in which the result needs to be stored
*/
__attribute__((always_inline)) INLINE static void
-rt_part_get_radiation_energy_density(const struct part *restrict p,
- float E[RT_NGROUPS]) {
+rt_part_get_comoving_radiation_energy_density(const struct part *restrict p,
+ float E[RT_NGROUPS]) {
for (int g = 0; g < RT_NGROUPS; g++) {
E[g] = p->rt_data.radiation[g].energy_density;
}
}
+/**
+ * @brief Get the physical radiation energy densities of a particle
+ *
+ * @param p Particle.
+ * @param E (return) Pointer to the array in which the result needs to be stored
+ */
+__attribute__((always_inline)) INLINE static void
+rt_part_get_physical_radiation_energy_density(const struct part *restrict p,
+ float E[RT_NGROUPS],
+ const struct cosmology *cosmo) {
+ for (int g = 0; g < RT_NGROUPS; g++) {
+ E[g] = cosmo->a3_inv * p->rt_data.radiation[g].energy_density;
+ }
+}
+
/**
* @brief Get a 4-element state vector U containing the radiation energy
* density and fluxes for a specific photon group.
diff --git a/src/rt/GEAR/rt_grackle_utils.h b/src/rt/GEAR/rt_grackle_utils.h
index 599221737365832fc1f8df665fe67757006a93e3..d5b281eabf3095769526798dc3298a855717b5c3 100644
--- a/src/rt/GEAR/rt_grackle_utils.h
+++ b/src/rt/GEAR/rt_grackle_utils.h
@@ -35,6 +35,21 @@
* @brief Utility and helper functions related to using grackle.
*/
+/**
+ * @brief Update grackle units during run
+ *
+ * @param grackle_units grackle units struct
+ * @param cosmo cosmology struct
+ *
+ * NOTE: In the current implementation, this function does nothing.
+ * However, there might be use-cases in the future (e.g. switching
+ * UV background on or off depending on redshift) that might be
+ * needed in the future, which can be implemented into this function.
+ */
+__attribute__((always_inline)) INLINE void update_grackle_units_cosmo(
+ code_units *grackle_units, const struct unit_system *us,
+ const struct cosmology *restrict cosmo) {}
+
/**
* @brief initialize grackle during rt_props_init
*
@@ -50,7 +65,8 @@ __attribute__((always_inline)) INLINE static void rt_init_grackle(
code_units *grackle_units, chemistry_data *grackle_chemistry_data,
chemistry_data_storage *grackle_chemistry_rates,
float hydrogen_mass_fraction, const int grackle_verb,
- const int case_B_recombination, const struct unit_system *us) {
+ const int case_B_recombination, const struct unit_system *us,
+ const struct cosmology *restrict cosmo) {
grackle_verbose = grackle_verb;
diff --git a/src/rt/GEAR/rt_properties.h b/src/rt/GEAR/rt_properties.h
index 49b0c9c60d50067f7927faa22bcc6df5f549db2d..6f9e95201bc8e145e11afc4c9f4a437814b6f205 100644
--- a/src/rt/GEAR/rt_properties.h
+++ b/src/rt/GEAR/rt_properties.h
@@ -452,7 +452,7 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
params, "GEARRT:case_B_recombination", /*default=*/1);
rt_init_grackle(&rtp->grackle_units, &rtp->grackle_chemistry_data,
&rtp->grackle_chemistry_rates, rtp->hydrogen_mass_fraction,
- rtp->grackle_verbose, rtp->case_B_recombination, us);
+ rtp->grackle_verbose, rtp->case_B_recombination, us, cosmo);
/* Pre-compute interaction rates/cross sections */
/* -------------------------------------------- */
@@ -465,6 +465,12 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
/* --------- */
}
+__attribute__((always_inline)) INLINE static void rt_props_update(
+ struct rt_props* rtp, const struct unit_system* us,
+ struct cosmology* cosmo) {
+ update_grackle_units_cosmo(&(rtp->grackle_units), us, cosmo);
+}
+
/**
* @brief Write an RT properties struct to the given FILE as a
* stream of bytes.
@@ -491,10 +497,11 @@ __attribute__((always_inline)) INLINE static void rt_struct_dump(
* @param stream the file stream
* @param phys_const The physical constants in the internal unit system.
* @param us The internal unit system.
+ * @param cosmo the #cosmology
*/
__attribute__((always_inline)) INLINE static void rt_struct_restore(
struct rt_props* props, FILE* stream, const struct phys_const* phys_const,
- const struct unit_system* us) {
+ const struct unit_system* us, const struct cosmology* restrict cosmo) {
restart_read_blocks((void*)props, sizeof(struct rt_props), 1, stream, NULL,
"RT properties struct");
@@ -502,7 +509,7 @@ __attribute__((always_inline)) INLINE static void rt_struct_restore(
rt_init_grackle(&props->grackle_units, &props->grackle_chemistry_data,
&props->grackle_chemistry_rates,
props->hydrogen_mass_fraction, props->grackle_verbose,
- props->case_B_recombination, us);
+ props->case_B_recombination, us, cosmo);
props->energy_weighted_cross_sections = NULL;
props->number_weighted_cross_sections = NULL;
diff --git a/src/rt/GEAR/rt_thermochemistry.h b/src/rt/GEAR/rt_thermochemistry.h
index 40aafd2b91f7e7f7c4a9a628d93c3b354bbcc352..a473f180393e601534c4225d8d18b72db62faf72 100644
--- a/src/rt/GEAR/rt_thermochemistry.h
+++ b/src/rt/GEAR/rt_thermochemistry.h
@@ -122,7 +122,6 @@ INLINE static void rt_do_thermochemistry(
grackle_field_data particle_grackle_data;
gr_float density = hydro_get_physical_density(p, cosmo);
-
/* In rare cases, unphysical solutions can arise with negative densities
* which won't be fixed in the hydro part until further down the dependency
* graph. Also, we can have vacuum, in which case we have nothing to do here.
@@ -130,15 +129,20 @@ INLINE static void rt_do_thermochemistry(
if (density <= 0.) return;
const float u_minimal = hydro_props->minimal_internal_energy;
- gr_float internal_energy =
- max(hydro_get_physical_internal_energy(p, xp, cosmo), u_minimal);
+
+ /* Physical internal energy */
+ gr_float internal_energy_phys =
+ hydro_get_physical_internal_energy(p, xp, cosmo);
+ gr_float internal_energy = max(internal_energy_phys, u_minimal);
+
const float u_old = internal_energy;
gr_float species_densities[6];
rt_tchem_get_species_densities(p, density, species_densities);
float radiation_energy_density[RT_NGROUPS];
- rt_part_get_radiation_energy_density(p, radiation_energy_density);
+ rt_part_get_physical_radiation_energy_density(p, radiation_energy_density,
+ cosmo);
gr_float iact_rates[5];
rt_get_interaction_rates_for_grackle(
@@ -152,9 +156,6 @@ INLINE static void rt_do_thermochemistry(
iact_rates);
/* solve chemistry */
- /* Note: `grackle_rates` is a global variable defined by grackle itself.
- * Using a manually allocd and initialized variable here fails with MPI
- * for some reason. */
if (local_solve_chemistry(
&rt_props->grackle_chemistry_data, &rt_props->grackle_chemistry_rates,
&rt_props->grackle_units, &particle_grackle_data, dt) == 0)
@@ -163,8 +164,9 @@ INLINE static void rt_do_thermochemistry(
/* copy updated grackle data to particle */
/* update particle internal energy. Grackle had access by reference
* to internal_energy */
- internal_energy = particle_grackle_data.internal_energy[0];
- const float u_new = max(internal_energy, u_minimal);
+ internal_energy_phys = particle_grackle_data.internal_energy[0];
+
+ const float u_new = max(internal_energy_phys, u_minimal);
/* Re-do thermochemistry? */
if ((rt_props->max_tchem_recursion > depth) &&
@@ -180,7 +182,7 @@ INLINE static void rt_do_thermochemistry(
}
/* If we're good, update the particle data from grackle results */
- hydro_set_internal_energy(p, u_new);
+ hydro_set_physical_internal_energy(p, xp, cosmo, u_new);
/* Update mass fractions */
const gr_float one_over_rho = 1. / density;
@@ -232,7 +234,6 @@ INLINE static void rt_do_thermochemistry(
p->rt_data.radiation[g].flux, E_old,
/*callloc=*/2);
}
-
/* Clean up after yourself. */
rt_clean_grackle_fields(&particle_grackle_data);
}
@@ -262,14 +263,18 @@ __attribute__((always_inline)) INLINE static float rt_tchem_get_tchem_time(
gr_float density = hydro_get_physical_density(p, cosmo);
const float u_minimal = hydro_props->minimal_internal_energy;
- gr_float internal_energy =
- max(hydro_get_physical_internal_energy(p, xp, cosmo), u_minimal);
+
+ /* Physical internal energy */
+ gr_float internal_energy_phys =
+ hydro_get_physical_internal_energy(p, xp, cosmo);
+ gr_float internal_energy = max(internal_energy_phys, u_minimal);
gr_float species_densities[6];
rt_tchem_get_species_densities(p, density, species_densities);
float radiation_energy_density[RT_NGROUPS];
- rt_part_get_radiation_energy_density(p, radiation_energy_density);
+ rt_part_get_physical_radiation_energy_density(p, radiation_energy_density,
+ cosmo);
gr_float iact_rates[5];
rt_get_interaction_rates_for_grackle(
@@ -282,15 +287,11 @@ __attribute__((always_inline)) INLINE static float rt_tchem_get_tchem_time(
iact_rates);
/* Compute 'cooling' time */
- /* Note: grackle_rates is a global variable defined by grackle itself.
- * Using a manually allocd and initialized variable here fails with MPI
- * for some reason. */
gr_float tchem_time;
if (local_calculate_cooling_time(
&rt_props->grackle_chemistry_data, &rt_props->grackle_chemistry_rates,
&rt_props->grackle_units, &particle_grackle_data, &tchem_time) == 0)
error("Error in calculate_cooling_time.");
-
/* Clean up after yourself. */
rt_clean_grackle_fields(&particle_grackle_data);
diff --git a/src/rt/SPHM1RT/rt_properties.h b/src/rt/SPHM1RT/rt_properties.h
index f2eb53d2030fb103296dd691fa3e207ac9fd2423..38dbde050a35a00a33eae3ef1f59637ca549ce49 100644
--- a/src/rt/SPHM1RT/rt_properties.h
+++ b/src/rt/SPHM1RT/rt_properties.h
@@ -402,6 +402,9 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
message("Radiative transfer initialized");
}
}
+__attribute__((always_inline)) INLINE static void rt_props_update(
+ struct rt_props* rtp, const struct unit_system* us,
+ struct cosmology* cosmo) {}
/**
* @brief Write an RT properties struct to the given FILE as a
@@ -428,7 +431,7 @@ __attribute__((always_inline)) INLINE static void rt_struct_dump(
*/
__attribute__((always_inline)) INLINE static void rt_struct_restore(
struct rt_props* props, FILE* stream, const struct phys_const* phys_const,
- const struct unit_system* us) {
+ const struct unit_system* us, const struct cosmology* restrict cosmo) {
restart_read_blocks((void*)props, sizeof(struct rt_props), 1, stream, NULL,
"RT properties struct");
diff --git a/src/rt/debug/rt_properties.h b/src/rt/debug/rt_properties.h
index 9fba0a49b274eb4499c31a6072032d92a15f45c8..d4a31a2c1b3c91537b3538f120bc1373ce41a0a8 100644
--- a/src/rt/debug/rt_properties.h
+++ b/src/rt/debug/rt_properties.h
@@ -92,6 +92,10 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
parser_get_param_int(params, "TimeIntegration:max_nr_rt_subcycles");
}
+__attribute__((always_inline)) INLINE static void rt_props_update(
+ struct rt_props* rtp, const struct unit_system* us,
+ struct cosmology* cosmo) {}
+
/**
* @brief Write an RT properties struct to the given FILE as a
* stream of bytes.
@@ -117,7 +121,7 @@ __attribute__((always_inline)) INLINE static void rt_struct_dump(
*/
__attribute__((always_inline)) INLINE static void rt_struct_restore(
struct rt_props* props, FILE* stream, const struct phys_const* phys_const,
- const struct unit_system* us) {
+ const struct unit_system* us, const struct cosmology* restrict cosmo) {
restart_read_blocks((void*)props, sizeof(struct rt_props), 1, stream, NULL,
"RT properties struct");
diff --git a/src/rt/none/rt_properties.h b/src/rt/none/rt_properties.h
index 16fd5fa5bf64987b7203e5437fa2da1ed1844aee..07e3cef313c2176b1d81f4a089472bd0b7c07f84 100644
--- a/src/rt/none/rt_properties.h
+++ b/src/rt/none/rt_properties.h
@@ -59,6 +59,10 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
const struct unit_system* us, struct swift_params* params,
struct cosmology* cosmo) {}
+__attribute__((always_inline)) INLINE static void rt_props_update(
+ struct rt_props* rtp, const struct unit_system* us,
+ struct cosmology* cosmo) {}
+
/**
* @brief Write an RT properties struct to the given FILE as a
* stream of bytes.
@@ -84,7 +88,7 @@ __attribute__((always_inline)) INLINE static void rt_struct_dump(
*/
__attribute__((always_inline)) INLINE static void rt_struct_restore(
struct rt_props* props, FILE* stream, const struct phys_const* phys_const,
- const struct unit_system* us) {
+ const struct unit_system* us, const struct cosmology* restrict cosmo) {
restart_read_blocks((void*)props, sizeof(struct rt_props), 1, stream, NULL,
"RT properties struct");