diff --git a/examples/RadiativeTransferTests/Advection_2D/makeIC.py b/examples/RadiativeTransferTests/Advection_2D/makeIC.py
index 923f9a56bf5964d9fa4dcea4854ce62e64a9fc04..e9742cbd6279129c696a6e5eb9230c29868870f4 100755
--- a/examples/RadiativeTransferTests/Advection_2D/makeIC.py
+++ b/examples/RadiativeTransferTests/Advection_2D/makeIC.py
@@ -40,7 +40,7 @@ from swiftsimio import Writer
# define unit system to use
unitsystem = unyt.unit_systems.cgs_unit_system
-# Box is 1 Mpc
+# define box size
boxsize = 1e10 * unitsystem["length"]
# number of photon groups
diff --git a/examples/RadiativeTransferTests/CollidingBeams_1D/README b/examples/RadiativeTransferTests/CollidingBeams_1D/README
new file mode 100644
index 0000000000000000000000000000000000000000..0aaee5bb32546bc83f47a0bbb1c9c3d9178fabea
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_1D/README
@@ -0,0 +1,13 @@
+Collide two photon packets with each other, and see what happens.
+Two photon groups contain packets using a top hat functions. One of
+them has the lower value at zero, the other is at nonzero. The third
+photon group packets are Gaussians.
+Ideally, they should just pass through each other. But that is not
+what moment based methods do with radiation.
+
+
+To run with GEAR-RT, compile SWIFT with
+
+ --with-rt=GEAR_2 --with-rt-riemann-solver=GLF --with-hydro-dimension=1 --with-hydro=gizmo-mfv \
+ --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
+
diff --git a/examples/RadiativeTransferTests/CollidingBeams_1D/makeIC.py b/examples/RadiativeTransferTests/CollidingBeams_1D/makeIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..ecfd7a24df2557f895a6115adffe9a6721f37c51
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_1D/makeIC.py
@@ -0,0 +1,192 @@
+#!/usr/bin/env python3
+
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2021 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+# 2022 Tsang Keung Chan (chantsangkeung@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/>.
+#
+##############################################################################
+
+
+# -----------------------------------------------------------
+# Add initial conditions for photon energies and fluxes
+# for 1D advection of photons.
+# First photon group: Top hat function with zero as the
+# baseline.
+# Second photon group: Top hat function with nonzero value
+# as the baseline.
+# Third photon group: Gaussian.
+# -----------------------------------------------------------
+
+import h5py
+import numpy as np
+import unyt
+from swiftsimio import Writer
+
+# define unit system to use
+unitsystem = unyt.unit_systems.cgs_unit_system
+
+# Box is 1 Mpc
+boxsize = 1e10 * unitsystem["length"]
+
+# number of photon groups
+nPhotonGroups = 3
+
+# number of particles in each dimension
+n_p = 1000
+
+# filename of ICs to be generated
+outputfilename = "collision_1D.hdf5"
+
+
+def initial_condition(x):
+ """
+ The initial conditions that will be advected
+
+ x: particle position. 3D unyt array
+
+ returns:
+ E: photon energy density for each photon group. List of scalars with size of nPhotonGroups
+ F: photon flux for each photon group. List with size of nPhotonGroups of numpy arrays of shape (3,)
+ """
+
+ # you can make the photon quantities unitless, the units will
+ # already have been written down in the writer.
+
+ E_list = []
+ F_list = []
+
+ # Group 1 Photons:
+ # -------------------
+
+ F = np.zeros(3, dtype=np.float64)
+
+ if x[0] < 0.33 * boxsize:
+ E = 1.0
+ F[0] = unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
+ elif x[0] < 0.66 * boxsize:
+ E = 0.0
+ F[0] = unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
+ else:
+ E = 1.0
+ F[0] = -unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
+
+ E_list.append(E)
+ F_list.append(F)
+
+ # Group 2 Photons:
+ # -------------------
+
+ F = np.zeros(3, dtype=np.float64)
+ if x[0] < 0.33 * boxsize:
+ E = 3.0
+ F[0] = unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
+ elif x[0] < 0.66 * boxsize:
+ E = 1.0
+ if x[0] < 0.5 * boxsize:
+ F[0] = unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
+ else:
+ F[0] = -unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
+ else:
+ E = 3.0
+ F[0] = -unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
+
+ E_list.append(E)
+ F_list.append(F)
+
+ # Group 3 Photons:
+ # -------------------
+ sigma = 0.05 * boxsize
+ if x[0] < 0.5 * boxsize:
+ mean = 0.33 / 2 * boxsize
+ sign = 1
+ else:
+ mean = (5.0 / 6.0) * boxsize
+ sign = -1
+ amplitude = 2.0
+
+ E = amplitude * np.exp(-(x[0] - mean) ** 2 / (2 * sigma ** 2))
+ F = np.zeros(3, dtype=np.float64)
+ F[0] = sign * unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
+
+ E_list.append(E)
+ F_list.append(F)
+
+ return E_list, F_list
+
+
+if __name__ == "__main__":
+
+ xp = unyt.unyt_array(np.zeros((n_p, 3), dtype=np.float64), boxsize.units)
+
+ dx = boxsize / n_p
+
+ for i in range(n_p):
+ xp[i, 0] = (i + 0.5) * dx
+
+ w = Writer(unyt.unit_systems.cgs_unit_system, boxsize, dimension=1)
+
+ w.gas.coordinates = xp
+ w.gas.velocities = np.zeros(xp.shape) * (unyt.cm / unyt.s)
+ w.gas.masses = np.ones(xp.shape[0], dtype=np.float64) * 1000 * unyt.g
+ w.gas.internal_energy = (
+ np.ones(xp.shape[0], dtype=np.float64) * (300.0 * unyt.kb * unyt.K) / unyt.g
+ )
+
+ # Generate initial guess for smoothing lengths based on MIPS
+ w.gas.generate_smoothing_lengths(boxsize=boxsize, dimension=1)
+
+ # If IDs are not present, this automatically generates
+ w.write(outputfilename)
+
+ # Now open file back up again and add photon groups
+ # you can make them unitless, the units have already been
+ # written down in the writer. In this case, it's in cgs.
+
+ F = h5py.File(outputfilename, "r+")
+ header = F["Header"]
+ nparts = header.attrs["NumPart_ThisFile"][0]
+ parts = F["/PartType0"]
+
+ for grp in range(nPhotonGroups):
+ dsetname = "PhotonEnergiesGroup{0:d}".format(grp + 1)
+ energydata = np.zeros((nparts), dtype=np.float32)
+ parts.create_dataset(dsetname, data=energydata)
+
+ dsetname = "PhotonFluxesGroup{0:d}".format(grp + 1)
+ # if dsetname not in parts.keys():
+ fluxdata = np.zeros((nparts, 3), dtype=np.float32)
+ parts.create_dataset(dsetname, data=fluxdata)
+
+ for p in range(nparts):
+ E, Flux = initial_condition(xp[p])
+ for g in range(nPhotonGroups):
+ Esetname = "PhotonEnergiesGroup{0:d}".format(g + 1)
+ parts[Esetname][p] = E[g]
+ Fsetname = "PhotonFluxesGroup{0:d}".format(g + 1)
+ parts[Fsetname][p] = Flux[g]
+
+ # from matplotlib import pyplot as plt
+ # plt.figure()
+ # for g in range(nPhotonGroups):
+ # Esetname = "PhotonEnergiesGroup{0:d}".format(g+1)
+ # plt.plot(xp[:,0], parts[Esetname], label="E "+str(g+1))
+ # # Fsetname = "PhotonFluxesGroup{0:d}".format(g+1)
+ # # plt.plot(xp[:,0], parts[Fsetname][:,0], label="F "+str(g+1))
+ # plt.legend()
+ # plt.show()
+
+ F.close()
diff --git a/examples/RadiativeTransferTests/CollidingBeams_1D/plotEnergies.py b/examples/RadiativeTransferTests/CollidingBeams_1D/plotEnergies.py
new file mode 100755
index 0000000000000000000000000000000000000000..4d5285467455e775e225a02a05378060070599c0
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_1D/plotEnergies.py
@@ -0,0 +1,159 @@
+#!/usr/bin/env python3
+
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Mladen Ivkovic (mladen.ivkovic@hotmail.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 total photon energies in each photon
+# group over the course of several snapshots
+# ----------------------------------------------
+
+import os
+import sys
+
+import matplotlib as mpl
+import numpy as np
+import swiftsimio
+from matplotlib import pyplot as plt
+
+# Parameters users should/may tweak
+snapshot_base = "output" # snapshot basename
+
+
+# -----------------------------------------------------------------------
+
+mpl.rcParams["text.usetex"] = True
+
+# Read in cmdline arg: Are we plotting only one snapshot, or all?
+plot_all = False # plot all snapshots
+try:
+ snapnr = int(sys.argv[1])
+except IndexError:
+ plot_all = True
+
+
+def get_snapshot_list(snapshot_basename="output"):
+ """
+ Find the snapshot(s) that are to be plotted
+ and return their names as list """
+
+ snaplist = []
+
+ if plot_all:
+ dirlist = os.listdir()
+ for f in dirlist:
+ if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+ snaplist.append(f)
+
+ snaplist = sorted(snaplist)
+
+ else:
+ fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5"
+ if not os.path.exists(fname):
+ print("Didn't find file", fname)
+ quit(1)
+ snaplist.append(fname)
+
+ return snaplist
+
+
+def plot_energies(snap_nrs, energy_sums):
+ """
+ Create the actual plot.
+ """
+
+ # Plot plot plot!
+ fig = plt.figure(figsize=(5.0, 5.4), dpi=200)
+ figname = "energy_budget.png"
+
+ ax1 = fig.add_subplot(1, 1, 1)
+
+ ngroups = energy_sums.shape[1]
+
+ for g in range(ngroups):
+ ax1.plot(snap_nrs, energy_sums[:, g], label=None)
+ ax1.scatter(snap_nrs, energy_sums[:, g], label="group {0:d}".format(g + 1))
+
+ ax1.set_xlabel("Snapshot")
+ ax1.set_ylabel(
+ r"Total energy [$" + energy_sums.units.latex_representation() + "$]",
+ usetex=True,
+ )
+
+ # add title
+ title = "Energy Budget"
+ ax1.set_title(title)
+ ax1.legend()
+
+ plt.tight_layout()
+ plt.savefig(figname)
+ plt.close()
+
+ return
+
+
+def get_photon_energies(snaplist):
+ """
+ Get total photon energy in each photon group for a list of
+ snapshots specified by `snaplist`
+
+ snaplist: list of snapshot filenames
+
+ returns:
+
+ snap_nrs : list of integers of snapshot numbers
+ energy_sums: np.array(shape=(len snaplist, ngroups)) of
+ total photon energies per group per snapshot
+ """
+
+ data = swiftsimio.load(snaplist[0])
+ meta = data.metadata
+ ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"])
+
+ energy_sums = np.zeros((len(snaplist), ngroups))
+ snap_nrs = np.zeros(len(snaplist), dtype=int)
+
+ for f, filename in enumerate(snaplist):
+
+ data = swiftsimio.load(filename)
+
+ for g in range(ngroups):
+ en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+ energy_sums[f, g] = en.sum()
+
+ nrstring = filename[len(snapshot_base) + 1 : -len(".hdf5")]
+ nr = int(nrstring)
+ snap_nrs[f] = nr
+
+ energy_sums = energy_sums * en.units
+
+ sortind = np.argsort(snap_nrs)
+ snap_nrs = snap_nrs[sortind]
+ energy_sums = energy_sums[sortind]
+
+ return snap_nrs, energy_sums
+
+
+if __name__ == "__main__":
+
+ snaplist = get_snapshot_list(snapshot_base)
+ snap_nrs, energy_sums = get_photon_energies(snaplist)
+
+ plot_energies(snap_nrs, energy_sums)
diff --git a/examples/RadiativeTransferTests/CollidingBeams_1D/plotSolution.py b/examples/RadiativeTransferTests/CollidingBeams_1D/plotSolution.py
new file mode 100755
index 0000000000000000000000000000000000000000..5e372d67976336086dac3579a524790242d238e8
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_1D/plotSolution.py
@@ -0,0 +1,264 @@
+#!/usr/bin/env python3
+
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Mladen Ivkovic (mladen.ivkovic@hotmail.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 photon data assuming a 1D problem
+# give snapshot number as cmdline arg to plot
+# single snapshot, otherwise this script plots
+# all snapshots available in the workdir
+# ----------------------------------------------
+
+import os
+import sys
+
+import matplotlib as mpl
+import numpy as np
+import swiftsimio
+from matplotlib import pyplot as plt
+
+# Parameters users should/may tweak
+snapshot_base = "output" # snapshot basename
+
+# properties for all scatterplots
+scatterplot_kwargs = {
+ "facecolor": "red",
+ "s": 4,
+ "alpha": 0.6,
+ "linewidth": 0.0,
+ "marker": ".",
+}
+
+# -----------------------------------------------------------------------
+
+mpl.rcParams["text.usetex"] = True
+
+# Read in cmdline arg: Are we plotting only one snapshot, or all?
+plot_all = False # plot all snapshots
+try:
+ snapnr = int(sys.argv[1])
+except IndexError:
+ plot_all = True
+
+
+def get_snapshot_list(snapshot_basename="output"):
+ """
+ Find the snapshot(s) that are to be plotted
+ and return their names as list
+ """
+
+ snaplist = []
+
+ if plot_all:
+ dirlist = os.listdir()
+ for f in dirlist:
+ if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+ snaplist.append(f)
+
+ snaplist = sorted(snaplist)
+
+ else:
+ fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5"
+ if not os.path.exists(fname):
+ print("Didn't find file", fname)
+ quit(1)
+ snaplist.append(fname)
+
+ return snaplist
+
+
+def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
+ """
+ Create the actual plot.
+
+ filename: file to work with
+ energy_boundaries: list of [E_min, E_max] for each photon group.
+ If none, limits are set automatically.
+ flux_boundaries: list of [F_min, F_max] for each photon group.
+ If none, limits are set automatically.
+ """
+
+ print("working on", filename)
+
+ # Read in data firt
+ data = swiftsimio.load(filename)
+ meta = data.metadata
+ scheme = str(meta.subgrid_scheme["RT Scheme"].decode("utf-8"))
+
+ ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"])
+
+ for g in range(ngroups):
+ # workaround to access named columns data with swiftsimio visualisaiton
+ new_attribute_str = "radiation_energy" + str(g + 1)
+ en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+ setattr(data.gas, new_attribute_str, en)
+
+ # prepare also the fluxes
+ for direction in ["X"]:
+ new_attribute_str = "radiation_flux" + str(g + 1) + direction
+ f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
+ setattr(data.gas, new_attribute_str, f)
+
+ part_positions = data.gas.coordinates[:, 0].copy()
+
+ # Plot plot plot!
+ fig = plt.figure(figsize=(5.0 * ngroups, 5.4), dpi=200)
+ figname = filename[:-5] + "-all-quantities.png"
+
+ for g in range(ngroups):
+
+ # plot energy
+ new_attribute_str = "radiation_energy" + str(g + 1)
+ photon_energy = getattr(data.gas, new_attribute_str)
+
+ ax = fig.add_subplot(2, ngroups, g + 1)
+
+ ax.scatter(
+ part_positions, photon_energy, **scatterplot_kwargs, label="simulation"
+ )
+ ax.legend()
+
+ ax.set_title("Group {0:2d}".format(g + 1))
+ if g == 0:
+ ax.set_ylabel(
+ "Energies [$" + photon_energy.units.latex_representation() + "$]"
+ )
+ ax.set_xlabel("x [$" + part_positions.units.latex_representation() + "$]")
+
+ if energy_boundaries is not None:
+
+ if abs(energy_boundaries[g][1]) > abs(energy_boundaries[g][0]):
+ fixed_min = energy_boundaries[g][0] - 0.1 * abs(energy_boundaries[g][1])
+ fixed_max = energy_boundaries[g][1] * 1.1
+ else:
+ fixed_min = energy_boundaries[g][0] * 1.1
+ fixed_max = energy_boundaries[g][1] + 0.1 * abs(energy_boundaries[g][0])
+ ax.set_ylim(fixed_min, fixed_max)
+
+ # plot flux X
+ new_attribute_str = "radiation_flux" + str(g + 1) + "X"
+ photon_flux = getattr(data.gas, new_attribute_str)
+ if scheme.startswith("GEAR M1closure"):
+ photon_flux = photon_flux.to("erg/cm**2/s")
+ elif scheme.startswith("SPH M1closure"):
+ photon_flux = photon_flux.to("erg*cm/s")
+ else:
+ print("Error: Unknown RT scheme " + scheme)
+ exit()
+
+ ax = fig.add_subplot(2, ngroups, g + 1 + ngroups)
+ ax.scatter(part_positions, photon_flux, **scatterplot_kwargs)
+
+ if g == 0:
+ ax.set_ylabel("Flux X [$" + photon_flux.units.latex_representation() + "$]")
+ ax.set_xlabel("x [$" + part_positions.units.latex_representation() + "$]")
+
+ if flux_boundaries is not None:
+
+ if abs(flux_boundaries[g][1]) > abs(flux_boundaries[g][0]):
+ fixed_min = flux_boundaries[g][0] - 0.1 * abs(flux_boundaries[g][1])
+ fixed_max = flux_boundaries[g][1] * 1.1
+ else:
+ fixed_min = flux_boundaries[g][0] * 1.1
+ fixed_max = flux_boundaries[g][1] + 0.1 * abs(flux_boundaries[g][0])
+
+ ax.set_ylim(fixed_min, fixed_max)
+
+ # add title
+ title = filename.replace("_", r"\_") # exception handle underscore for latex
+ if meta.cosmology is not None:
+ title += ", $z$ = {0:.2e}".format(meta.z)
+ title += ", $t$ = {0:.2e}".format(meta.time)
+ fig.suptitle(title)
+
+ plt.tight_layout()
+ plt.savefig(figname)
+ plt.close()
+
+ return
+
+
+def get_minmax_vals(snaplist):
+ """
+ Find minimal and maximal values for energy and flux,
+ so you can fix axes limits over all snapshots
+
+ snaplist: list of snapshot filenames
+
+ returns:
+
+ energy_boundaries: list of [E_min, E_max] for each photon group
+ flux_boundaries: list of [Fx_min, Fy_max] for each photon group
+ """
+
+ emins = []
+ emaxs = []
+ fmins = []
+ fmaxs = []
+
+ for filename in snaplist:
+
+ data = swiftsimio.load(filename)
+ meta = data.metadata
+
+ ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"])
+ emin_group = []
+ emax_group = []
+ fluxmin_group = []
+ fluxmax_group = []
+
+ for g in range(ngroups):
+ en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+ emin_group.append(en.min())
+ emax_group.append(en.max())
+
+ for direction in ["X"]:
+ f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
+ fluxmin_group.append(f.min())
+ fluxmax_group.append(f.max())
+
+ emins.append(emin_group)
+ emaxs.append(emax_group)
+ fmins.append(fluxmin_group)
+ fmaxs.append(fluxmax_group)
+
+ energy_boundaries = []
+ flux_boundaries = []
+ for g in range(ngroups):
+ emin = min([emins[f][g] for f in range(len(snaplist))])
+ emax = max([emaxs[f][g] for f in range(len(snaplist))])
+ energy_boundaries.append([emin, emax])
+ fmin = min([fmins[f][g] for f in range(len(snaplist))])
+ fmax = max([fmaxs[f][g] for f in range(len(snaplist))])
+ flux_boundaries.append([fmin, fmax])
+
+ return energy_boundaries, flux_boundaries
+
+
+if __name__ == "__main__":
+
+ snaplist = get_snapshot_list(snapshot_base)
+ energy_boundaries, flux_boundaries = get_minmax_vals(snaplist)
+
+ for f in snaplist:
+ plot_photons(
+ f, energy_boundaries=energy_boundaries, flux_boundaries=flux_boundaries
+ )
diff --git a/examples/RadiativeTransferTests/CollidingBeams_1D/rt_collision1D.yml b/examples/RadiativeTransferTests/CollidingBeams_1D/rt_collision1D.yml
new file mode 100644
index 0000000000000000000000000000000000000000..25f24418dada5071893d5b8bc1679dfd77c4b614
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_1D/rt_collision1D.yml
@@ -0,0 +1,62 @@
+MetaData:
+ run_name: RT_collision-1D
+
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.
+ UnitLength_in_cgs: 1.
+ UnitVelocity_in_cgs: 1.
+ UnitCurrent_in_cgs: 1.
+ UnitTemp_in_cgs: 1.
+
+# 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: 4.e-1 # 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: 1.e-02 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: output # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 2.e-2
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ time_first: 0.
+ delta_time: 4.e-2 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ CFL_condition: 0.6 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 10. # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./collision_1D.hdf5 # The file to read
+ periodic: 1 # periodic ICs
+
+Restarts:
+ delta_hours: 72 # (Optional) decimal hours between dumps of restart files.
+
+GEARRT:
+ f_reduce_c: 1. # reduce the speed of light for the RT solver by multiplying c with this factor
+ f_limit_cooling_time: 0.0 # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+ CFL_condition: 0.8 # CFL condition for RT, independent of hydro
+ photon_groups_Hz: [0., 1., 2.] # 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 (H + H+) mass fraction in the metal-free portion of the gas
+ 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.
+ skip_thermochemistry: 1 # ignore thermochemistry.
+ set_initial_ionization_mass_fractions: 1 # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+ mass_fraction_HI: 0.76 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+ mass_fraction_HII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+ mass_fraction_HeI: 0.24 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+ mass_fraction_HeII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+ mass_fraction_HeIII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+
diff --git a/examples/RadiativeTransferTests/CollidingBeams_1D/run.sh b/examples/RadiativeTransferTests/CollidingBeams_1D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..1d121bcc7eeefbb6e257057f78cae61e1981cf58
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_1D/run.sh
@@ -0,0 +1,24 @@
+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -f collision_1D.hdf5 ]; then
+ echo "Generating ICs"
+ python3 makeIC.py
+fi
+
+# Run SWIFT with RT
+../../../swift \
+ --hydro \
+ --threads=1 \
+ --verbose=0 \
+ --radiation \
+ --stars \
+ --feedback \
+ --external-gravity \
+ ./rt_collision1D.yml 2>&1 | tee output.log
+
+python3 ./plotEnergies.py
+python3 ./plotSolution.py
diff --git a/examples/RadiativeTransferTests/CollidingBeams_2D/README b/examples/RadiativeTransferTests/CollidingBeams_2D/README
new file mode 100644
index 0000000000000000000000000000000000000000..7bc95ce3b54b0281536b7c1fea3bbbfce267e4e9
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_2D/README
@@ -0,0 +1,10 @@
+Collide packets of radiation against each other, and see how the M1
+closure treats them.
+
+Collide them horizontally, vertically, diagonally, and from the side.
+
+To run with GEAR-RT, configure SWIFT using
+
+ --with-rt=GEAR_4 --with-rt-riemann-solver=GLF --with-hydro-dimension=2 --with-hydro=gizmo-mfv \
+ --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
+
diff --git a/examples/RadiativeTransferTests/CollidingBeams_2D/getGlass.sh b/examples/RadiativeTransferTests/CollidingBeams_2D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ae3c977064f5e7a408aa249c5fd9089b3c52ecb1
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_2D/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_128.hdf5
diff --git a/examples/RadiativeTransferTests/CollidingBeams_2D/makeIC.py b/examples/RadiativeTransferTests/CollidingBeams_2D/makeIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..be2d50a7bb2a6503195a0d6c700f3ee6615091de
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_2D/makeIC.py
@@ -0,0 +1,302 @@
+#!/usr/bin/env python3
+
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Mladen Ivkovic (mladen.ivkovic@hotmail.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/>.
+#
+##############################################################################
+
+
+# -------------------------------------------------------------
+# Add initial conditions for photon energies and fluxes
+# for 2D advection of photons.
+# First photon group: Top hat function with zero as the
+# baseline, advects in x direction
+# Second photon group: Top hat function with nonzero value
+# as the baseline, advcts in y direction.
+# Third photon group: Gaussian advecting diagonally
+# Fourth photon group: Circle moving radially from the center
+# -------------------------------------------------------------
+
+import h5py
+import numpy as np
+import unyt
+from swiftsimio import Writer
+
+# define unit system to use
+unitsystem = unyt.unit_systems.cgs_unit_system
+
+# define box size
+boxsize = 1e10 * unitsystem["length"]
+
+# number of photon groups
+nPhotonGroups = 4
+
+# filename of ICs to be generated
+outputfilename = "collision_2D.hdf5"
+
+
+def get_radiation_IC(pos, numPart, ngroups):
+
+ Elist = []
+ Flist = []
+
+ x = pos[:, 0]
+ y = pos[:, 1]
+ c = unyt.c.to(unitsystem["length"] / unitsystem["time"])
+
+ for g in range(ngroups):
+
+ E = np.zeros(numPart)
+ F = np.zeros((numPart, 3))
+
+ # prepare particle array masks
+ flux_sign = np.ones(numPart, dtype=float)
+
+ if g == 0:
+ # horizontal beams
+ # -----------------------
+
+ left = np.logical_and(x > 0.0, x < 1.0 / 3)
+ right = np.logical_and(x > 2.0 / 3, x < 1.0)
+ xcriterion = np.logical_or(left, right)
+ ycriterion = np.logical_and(y > 0.4, y < 0.6)
+ is_max = np.logical_and(xcriterion, ycriterion)
+
+ flux_sign[right] = -1.0
+
+ Emax = 1.0
+ E[is_max] = Emax
+ F[is_max, 0] = Emax * c * flux_sign[is_max]
+
+ if g == 1:
+ # vertical beams, nonzero energy everywhere
+ # -------------------------------------------
+
+ top = np.logical_and(y > 2.0 / 3.0, y < 1.0)
+ bottom = np.logical_and(y > 0.0, y < 1.0 / 3.0)
+
+ xcriterion = np.logical_and(x > 0.4, x < 0.6)
+ ycriterion = np.logical_or(top, bottom)
+ is_max = np.logical_and(xcriterion, ycriterion)
+
+ flux_sign[y > 0.5] = -1.0
+
+ Emax = 2.0
+ Emin = 1.0
+ E[:] = Emin
+ E[is_max] = Emax
+ F[:, 1] = E * c * flux_sign
+
+ if g == 2:
+ # diagonal beams
+ # -------------------------------------------
+
+ width = 0.1
+ l = np.sqrt(2) / 3.0
+ sin = np.sin(np.pi / 4)
+ xprime = l * sin
+ x0 = xprime - 0.5 * width * sin
+ x1 = xprime + 0.5 * width * sin
+ x2 = 1 - xprime - 0.5 * width * sin
+ x3 = 1 - xprime + 0.5 * width * sin
+
+ def upper_line(x):
+ return x + 0.5 * width
+
+ def lower_line(x):
+ return x - 0.5 * width
+
+ def descending_left(x, xprime):
+ return -(x - xprime) + xprime
+
+ def descending_right(x, xprime):
+ return -(x - 1.0 + xprime) + 1.0 - xprime
+
+ first = x < x0
+ lower_first = np.logical_and(y < upper_line(x), y > lower_line(x))
+ firstcond = np.logical_and(first, lower_first)
+
+ second = np.logical_and(x > x0, x < x1)
+ lower_second = np.logical_and(
+ y > lower_line(x), y < descending_left(x, xprime)
+ )
+ secondcond = np.logical_and(second, lower_second)
+
+ third = np.logical_and(x > x2, x < x3)
+ upper_third = np.logical_and(
+ y < upper_line(x), y > descending_right(x, xprime)
+ )
+ thirdcond = np.logical_and(third, upper_third)
+
+ fourth = np.logical_and(x > x3, x < 1.0)
+ upper_fourth = np.logical_and(y < upper_line(x), y > lower_line(x))
+ fourthcond = np.logical_and(fourth, upper_fourth)
+
+ Emax = 1.0
+ E[firstcond] = Emax
+ E[secondcond] = Emax
+ E[thirdcond] = Emax
+ E[fourthcond] = Emax
+
+ flux_sign[thirdcond] = -1.0
+ flux_sign[fourthcond] = -1.0
+
+ F[:, 0] = E * c * flux_sign / np.sqrt(2)
+ F[:, 1] = E * c * flux_sign / np.sqrt(2)
+
+ if g == 3:
+ # diagonal beams that meed in the middle
+ # -------------------------------------------
+
+ width = 0.1
+ l = np.sqrt(2) / 3.0
+ sin = np.sin(np.pi / 4)
+ xprime = l * sin
+ x0 = xprime - 0.5 * width * sin
+ x1 = xprime + 0.5 * width * sin
+ x2 = 1 - xprime - 0.5 * width * sin
+ x3 = 1 - xprime + 0.5 * width * sin
+
+ def upper_line(x):
+ return x + 0.5 * width
+
+ def lower_line(x):
+ return x - 0.5 * width
+
+ def descending_left(x, xprime):
+ return -(x - xprime) + xprime
+
+ def descending_lower(x, xprime):
+ return -(x - 1 + xprime) + xprime - 0.5 * width
+
+ def descending_upper(x, xprime):
+ return -(x - 1.0 + xprime) + xprime + 0.5 * width
+
+ def ascending_right(x, xprime):
+ return (x - 1 + xprime) + xprime
+
+ first = x < x0
+ lower_first = np.logical_and(y < upper_line(x), y > lower_line(x))
+ firstcond = np.logical_and(first, lower_first)
+
+ second = np.logical_and(x > x0, x < x1)
+ lower_second = np.logical_and(
+ y > lower_line(x), y < descending_left(x, xprime)
+ )
+ secondcond = np.logical_and(second, lower_second)
+
+ third = np.logical_and(x > x2, x < x3)
+ upper_third = np.logical_and(
+ y < ascending_right(x, xprime), y > descending_lower(x, xprime)
+ )
+ thirdcond = np.logical_and(third, upper_third)
+
+ fourth = np.logical_and(x > x3, x < 1.0)
+ upper_fourth = np.logical_and(
+ y < descending_upper(x, xprime), y > descending_lower(x, xprime)
+ )
+ fourthcond = np.logical_and(fourth, upper_fourth)
+
+ Emax = 1.0
+ E[firstcond] = Emax
+ E[secondcond] = Emax
+ E[thirdcond] = Emax
+ E[fourthcond] = Emax
+
+ flux_sign[thirdcond] = -1.0
+ flux_sign[fourthcond] = -1.0
+
+ F[:, 0] = E * c * flux_sign / np.sqrt(2)
+ F[:, 1] = E * c / np.sqrt(2)
+
+ # histE, _, _ = np.histogram2d(pos[:,0], pos[:,1], weights=E, bins=50)
+ # histFx, _, _ = np.histogram2d(pos[:,0], pos[:,1], weights=F[:,0], bins=50)
+ # histFy, _, _ = np.histogram2d(pos[:,0], pos[:,1], weights=F[:,1], bins=50)
+ # from matplotlib import pyplot as plt
+ #
+ # fig = plt.figure()
+ # ax1 = fig.add_subplot(1, 3, 1)
+ # ax1.imshow(histE.T, origin="lower")
+ # ax2 = fig.add_subplot(1, 3, 2)
+ # ax2.imshow(histFx.T, origin="lower")
+ # ax3 = fig.add_subplot(1, 3, 3)
+ # ax3.imshow(histFy.T, origin="lower")
+ # plt.show()
+
+ Elist.append(E)
+ Flist.append(F)
+
+ return Elist, Flist
+
+
+if __name__ == "__main__":
+ glass = h5py.File("glassPlane_128.hdf5", "r")
+
+ # Read particle positions and h from the glass
+ pos = glass["/PartType0/Coordinates"][:, :]
+ h = glass["/PartType0/SmoothingLength"][:]
+ glass.close()
+
+ numPart = np.size(h)
+
+ # get radiation IC while particle coordinates are still [0, 1)
+ Elist, Flist = get_radiation_IC(pos, numPart, nPhotonGroups)
+
+ pos *= boxsize
+ h *= boxsize
+
+ w = Writer(unyt.unit_systems.cgs_unit_system, boxsize, dimension=2)
+
+ w.gas.coordinates = pos
+ w.gas.velocities = np.zeros((numPart, 3)) * (unyt.cm / unyt.s)
+ w.gas.masses = np.ones(numPart, dtype=np.float64) * 1000 * unyt.g
+ w.gas.internal_energy = (
+ np.ones(numPart, dtype=np.float64) * (300.0 * unyt.kb * unyt.K) / unyt.g
+ )
+
+ # Generate initial guess for smoothing lengths based on MIPS
+ w.gas.smoothing_length = h
+
+ # If IDs are not present, this automatically generates
+ w.write(outputfilename)
+
+ # Now open file back up again and add photon groups
+ # you can make them unitless, the units have already been
+ # written down in the writer. In this case, it's in cgs.
+
+ F = h5py.File(outputfilename, "r+")
+ header = F["Header"]
+ nparts = header.attrs["NumPart_ThisFile"][0]
+ parts = F["/PartType0"]
+
+ for grp in range(nPhotonGroups):
+ dsetname = "PhotonEnergiesGroup{0:d}".format(grp + 1)
+ energydata = np.zeros(nparts, dtype=np.float32)
+ parts.create_dataset(dsetname, data=energydata)
+
+ dsetname = "PhotonFluxesGroup{0:d}".format(grp + 1)
+ # if dsetname not in parts.keys():
+ fluxdata = np.zeros((nparts, 3), dtype=np.float32)
+ parts.create_dataset(dsetname, data=fluxdata)
+
+ for g in range(nPhotonGroups):
+ Esetname = "PhotonEnergiesGroup{0:d}".format(g + 1)
+ parts[Esetname][:] = Elist[g][:]
+ Fsetname = "PhotonFluxesGroup{0:d}".format(g + 1)
+ parts[Fsetname][:, :] = Flist[g][:, :]
+
+ F.close()
diff --git a/examples/RadiativeTransferTests/CollidingBeams_2D/plotEnergies.py b/examples/RadiativeTransferTests/CollidingBeams_2D/plotEnergies.py
new file mode 100755
index 0000000000000000000000000000000000000000..4d5285467455e775e225a02a05378060070599c0
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_2D/plotEnergies.py
@@ -0,0 +1,159 @@
+#!/usr/bin/env python3
+
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Mladen Ivkovic (mladen.ivkovic@hotmail.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 total photon energies in each photon
+# group over the course of several snapshots
+# ----------------------------------------------
+
+import os
+import sys
+
+import matplotlib as mpl
+import numpy as np
+import swiftsimio
+from matplotlib import pyplot as plt
+
+# Parameters users should/may tweak
+snapshot_base = "output" # snapshot basename
+
+
+# -----------------------------------------------------------------------
+
+mpl.rcParams["text.usetex"] = True
+
+# Read in cmdline arg: Are we plotting only one snapshot, or all?
+plot_all = False # plot all snapshots
+try:
+ snapnr = int(sys.argv[1])
+except IndexError:
+ plot_all = True
+
+
+def get_snapshot_list(snapshot_basename="output"):
+ """
+ Find the snapshot(s) that are to be plotted
+ and return their names as list """
+
+ snaplist = []
+
+ if plot_all:
+ dirlist = os.listdir()
+ for f in dirlist:
+ if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+ snaplist.append(f)
+
+ snaplist = sorted(snaplist)
+
+ else:
+ fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5"
+ if not os.path.exists(fname):
+ print("Didn't find file", fname)
+ quit(1)
+ snaplist.append(fname)
+
+ return snaplist
+
+
+def plot_energies(snap_nrs, energy_sums):
+ """
+ Create the actual plot.
+ """
+
+ # Plot plot plot!
+ fig = plt.figure(figsize=(5.0, 5.4), dpi=200)
+ figname = "energy_budget.png"
+
+ ax1 = fig.add_subplot(1, 1, 1)
+
+ ngroups = energy_sums.shape[1]
+
+ for g in range(ngroups):
+ ax1.plot(snap_nrs, energy_sums[:, g], label=None)
+ ax1.scatter(snap_nrs, energy_sums[:, g], label="group {0:d}".format(g + 1))
+
+ ax1.set_xlabel("Snapshot")
+ ax1.set_ylabel(
+ r"Total energy [$" + energy_sums.units.latex_representation() + "$]",
+ usetex=True,
+ )
+
+ # add title
+ title = "Energy Budget"
+ ax1.set_title(title)
+ ax1.legend()
+
+ plt.tight_layout()
+ plt.savefig(figname)
+ plt.close()
+
+ return
+
+
+def get_photon_energies(snaplist):
+ """
+ Get total photon energy in each photon group for a list of
+ snapshots specified by `snaplist`
+
+ snaplist: list of snapshot filenames
+
+ returns:
+
+ snap_nrs : list of integers of snapshot numbers
+ energy_sums: np.array(shape=(len snaplist, ngroups)) of
+ total photon energies per group per snapshot
+ """
+
+ data = swiftsimio.load(snaplist[0])
+ meta = data.metadata
+ ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"])
+
+ energy_sums = np.zeros((len(snaplist), ngroups))
+ snap_nrs = np.zeros(len(snaplist), dtype=int)
+
+ for f, filename in enumerate(snaplist):
+
+ data = swiftsimio.load(filename)
+
+ for g in range(ngroups):
+ en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+ energy_sums[f, g] = en.sum()
+
+ nrstring = filename[len(snapshot_base) + 1 : -len(".hdf5")]
+ nr = int(nrstring)
+ snap_nrs[f] = nr
+
+ energy_sums = energy_sums * en.units
+
+ sortind = np.argsort(snap_nrs)
+ snap_nrs = snap_nrs[sortind]
+ energy_sums = energy_sums[sortind]
+
+ return snap_nrs, energy_sums
+
+
+if __name__ == "__main__":
+
+ snaplist = get_snapshot_list(snapshot_base)
+ snap_nrs, energy_sums = get_photon_energies(snaplist)
+
+ plot_energies(snap_nrs, energy_sums)
diff --git a/examples/RadiativeTransferTests/CollidingBeams_2D/plotSolution.py b/examples/RadiativeTransferTests/CollidingBeams_2D/plotSolution.py
new file mode 100755
index 0000000000000000000000000000000000000000..8f100e5fae99c15fb18d63b3f572db6f045d802f
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_2D/plotSolution.py
@@ -0,0 +1,310 @@
+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2021 Mladen Ivkovic (mladen.ivkovic@hotmail.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 photon data for 2D problems
+# give snapshot number as cmdline arg to plot
+# single snapshot, otherwise this script plots
+# all snapshots available in the workdir
+# ----------------------------------------------------
+
+import gc
+import os
+import sys
+
+import matplotlib as mpl
+import swiftsimio
+from matplotlib import pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+# Parameters users should/may tweak
+plot_all_data = False # plot all groups and all photon quantities
+snapshot_base = "output" # snapshot basename
+fancy = True # fancy up the plots a bit?
+
+# parameters for imshow plots
+imshow_kwargs = {"origin": "lower", "cmap": "viridis"}
+
+
+projection_kwargs = {"resolution": 1024, "parallel": True}
+# -----------------------------------------------------------------------
+
+
+# Read in cmdline arg: Are we plotting only one snapshot, or all?
+plot_all = False
+try:
+ snapnr = int(sys.argv[1])
+except IndexError:
+ plot_all = True
+
+mpl.rcParams["text.usetex"] = True
+
+
+def get_snapshot_list(snapshot_basename="output"):
+ """
+ Find the snapshot(s) that are to be plotted
+ and return their names as list
+ """
+
+ snaplist = []
+
+ if plot_all:
+ dirlist = os.listdir()
+ for f in dirlist:
+ if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+ snaplist.append(f)
+
+ snaplist = sorted(snaplist)
+
+ else:
+ fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5"
+ if not os.path.exists(fname):
+ print("Didn't find file", fname)
+ quit(1)
+ snaplist.append(fname)
+
+ return snaplist
+
+
+def set_colorbar(ax, im):
+ divider = make_axes_locatable(ax)
+ cax = divider.append_axes("right", size="5%", pad=0.05)
+ plt.colorbar(im, cax=cax)
+ return
+
+
+def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
+ """
+ Create the actual plot.
+
+ filename: file to work with
+ energy_boundaries: list of [E_min, E_max] for each photon group.
+ If none, limits are set automatically.
+ flux_boundaries: list of [F_min, F_max] for each photon group.
+ If none, limits are set automatically.
+ """
+
+ print("working on", filename)
+
+ # Read in data first
+ data = swiftsimio.load(filename)
+ meta = data.metadata
+
+ ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"])
+ xlabel_units_str = meta.boxsize.units.latex_representation()
+
+ global imshow_kwargs
+ imshow_kwargs["extent"] = [0, meta.boxsize[0].v, 0, meta.boxsize[1].v]
+
+ for g in range(ngroups):
+ # workaround to access named columns data with swiftsimio visualisaiton
+ # add mass weights to remove surface density dependence in images
+ new_attribute_str = "mass_weighted_radiation_energy" + str(g + 1)
+ en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+ en = en * data.gas.masses
+ setattr(data.gas, new_attribute_str, en)
+
+ if plot_all_data:
+ # prepare also the fluxes
+ # for direction in ["X", "Y", "Z"]:
+ for direction in ["X", "Y"]:
+ new_attribute_str = (
+ "mass_weighted_radiation_flux" + str(g + 1) + direction
+ )
+ f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
+ f *= data.gas.masses
+ setattr(data.gas, new_attribute_str, f)
+
+ # get mass surface density projection that we'll use to remove density dependence in image
+ mass_map = swiftsimio.visualisation.projection.project_gas(
+ data, project="masses", **projection_kwargs
+ )
+
+ if plot_all_data:
+ fig = plt.figure(figsize=(5 * 3, 5.05 * ngroups), dpi=200)
+ figname = filename[:-5] + "-all-quantities.png"
+
+ for g in range(ngroups):
+
+ # get energy projection
+ new_attribute_str = "mass_weighted_radiation_energy" + str(g + 1)
+ photon_map = swiftsimio.visualisation.projection.project_gas(
+ data, project=new_attribute_str, **projection_kwargs
+ )
+ photon_map = photon_map / mass_map
+
+ ax = fig.add_subplot(ngroups, 3, g * 3 + 1)
+ if energy_boundaries is not None:
+ imshow_kwargs["vmin"] = energy_boundaries[g][0]
+ imshow_kwargs["vmax"] = energy_boundaries[g][1]
+ im = ax.imshow(photon_map.T, **imshow_kwargs)
+ set_colorbar(ax, im)
+ ax.set_ylabel("Group {0:2d}".format(g + 1))
+ ax.set_xlabel("x [$" + xlabel_units_str + "$]")
+ if g == 0:
+ ax.set_title("Energies")
+
+ # get flux X projection
+ new_attribute_str = "mass_weighted_radiation_flux" + str(g + 1) + "X"
+ photon_map = swiftsimio.visualisation.projection.project_gas(
+ data, project=new_attribute_str, **projection_kwargs
+ )
+ photon_map = photon_map / mass_map
+
+ ax = fig.add_subplot(ngroups, 3, g * 3 + 2)
+ if flux_boundaries is not None:
+ imshow_kwargs["vmin"] = flux_boundaries[g][0]
+ imshow_kwargs["vmax"] = flux_boundaries[g][1]
+ im = ax.imshow(photon_map.T, **imshow_kwargs)
+ set_colorbar(ax, im)
+ ax.set_xlabel("x [$" + xlabel_units_str + "$]")
+ ax.set_ylabel("y [$" + xlabel_units_str + "$]")
+ if g == 0:
+ ax.set_title("Flux X")
+
+ # get flux Y projection
+ new_attribute_str = "mass_weighted_radiation_flux" + str(g + 1) + "Y"
+ photon_map = swiftsimio.visualisation.projection.project_gas(
+ data, project=new_attribute_str, **projection_kwargs
+ )
+ photon_map = photon_map / mass_map
+
+ ax = fig.add_subplot(ngroups, 3, g * 3 + 3)
+ im = ax.imshow(photon_map.T, **imshow_kwargs)
+ set_colorbar(ax, im)
+ ax.set_xlabel("x [$" + xlabel_units_str + "$]")
+ ax.set_ylabel("y [$" + xlabel_units_str + "$]")
+ if g == 0:
+ ax.set_title("Flux Y")
+
+ else: # plot just energies
+
+ fig = plt.figure(figsize=(5 * ngroups, 5), dpi=200)
+ figname = filename[:-5] + ".png"
+
+ for g in range(ngroups):
+
+ # get projection
+ new_attribute_str = "mass_weighted_radiation_energy" + str(g + 1)
+ photon_map = swiftsimio.visualisation.projection.project_gas(
+ data, project=new_attribute_str, **projection_kwargs
+ )
+ photon_map = photon_map / mass_map
+
+ ax = fig.add_subplot(1, ngroups, g + 1)
+ if energy_boundaries is not None:
+ imshow_kwargs["vmin"] = energy_boundaries[g][0]
+ imshow_kwargs["vmax"] = energy_boundaries[g][1]
+ im = ax.imshow(photon_map.T, **imshow_kwargs)
+ set_colorbar(ax, im)
+ ax.set_title("Group {0:2d}".format(g + 1))
+ if g == 0:
+ ax.set_ylabel("Energies")
+
+ # Add title
+ title = filename.replace("_", r"\_") # exception handle underscore for latex
+ if meta.cosmology is not None:
+ title += ", $z$ = {0:.2e}".format(meta.z)
+ title += ", $t$ = {0:.2e}".format(meta.time)
+ fig.suptitle(title)
+
+ plt.tight_layout()
+ plt.savefig(figname)
+ plt.close()
+ gc.collect()
+
+ return
+
+
+def get_minmax_vals(snaplist):
+ """
+ Find minimal and maximal values for energy and flux,
+ so you can fix axes limits over all snapshots
+
+ snaplist: list of snapshot filenames
+
+ returns:
+
+ energy_boundaries: list of [E_min, E_max] for each photon group
+ flux_boundaries: list of [Fx_min, Fy_max] for each photon group
+ """
+
+ emins = []
+ emaxs = []
+ fmins = []
+ fmaxs = []
+
+ for filename in snaplist:
+
+ data = swiftsimio.load(filename)
+ meta = data.metadata
+
+ ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"])
+ emin_group = []
+ emax_group = []
+ fluxmin_group = []
+ fluxmax_group = []
+
+ for g in range(ngroups):
+ en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+ emin_group.append(en.min())
+ emax_group.append(en.max())
+
+ dirmin = []
+ dirmax = []
+ for direction in ["X", "Y"]:
+ f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
+ dirmin.append(f.min())
+ dirmax.append(f.max())
+ fluxmin_group.append(min(dirmin))
+ fluxmax_group.append(max(dirmax))
+
+ emins.append(emin_group)
+ emaxs.append(emax_group)
+ fmins.append(fluxmin_group)
+ fmaxs.append(fluxmax_group)
+
+ energy_boundaries = []
+ flux_boundaries = []
+ for g in range(ngroups):
+ emin = min([emins[f][g] for f in range(len(snaplist))])
+ emax = max([emaxs[f][g] for f in range(len(snaplist))])
+ energy_boundaries.append([emin, emax])
+ fmin = min([fmins[f][g] for f in range(len(snaplist))])
+ fmax = max([fmaxs[f][g] for f in range(len(snaplist))])
+ flux_boundaries.append([fmin, fmax])
+
+ return energy_boundaries, flux_boundaries
+
+
+if __name__ == "__main__":
+
+ snaplist = get_snapshot_list(snapshot_base)
+ if fancy:
+ energy_boundaries, flux_boundaries = get_minmax_vals(snaplist)
+ else:
+ energy_boundaries = None
+ flux_boundaries = None
+
+ for f in snaplist:
+ plot_photons(
+ f, energy_boundaries=energy_boundaries, flux_boundaries=flux_boundaries
+ )
diff --git a/examples/RadiativeTransferTests/CollidingBeams_2D/rt_collision2D.yml b/examples/RadiativeTransferTests/CollidingBeams_2D/rt_collision2D.yml
new file mode 100644
index 0000000000000000000000000000000000000000..44daff0259e0bd8165b8d18beb13afe459bb040c
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_2D/rt_collision2D.yml
@@ -0,0 +1,60 @@
+MetaData:
+ run_name: RT_collision-2D
+
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.
+ UnitLength_in_cgs: 1.
+ UnitVelocity_in_cgs: 1.
+ UnitCurrent_in_cgs: 1.
+ UnitTemp_in_cgs: 1.
+
+# 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: 3.4e-1 # The end time of the simulation (in internal units).
+ dt_min: 1.e-08 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1.e-02 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: output # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.017
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ time_first: 0.
+ delta_time: 0.034 # 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: 10. # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./collision_2D.hdf5 # The file to read
+ periodic: 1 # periodic ICs
+
+GEARRT:
+ f_reduce_c: 1.
+ # reduce the speed of light for the RT solver by multiplying c with this factor
+ f_limit_cooling_time: 0.0 # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+ CFL_condition: 0.8 # CFL condition for RT, independent of hydro
+ photon_groups_Hz: [1., 2., 3., 4.] # 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., 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 (H + H+) mass fraction in the metal-free portion of the gas
+ 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.
+ set_initial_ionization_mass_fractions: 1 # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+ mass_fraction_HI: 0.76 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+ mass_fraction_HII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+ mass_fraction_HeI: 0.24 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+ mass_fraction_HeII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+ mass_fraction_HeIII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+ skip_thermochemistry: 1 # Skip thermochemsitry.
+
diff --git a/examples/RadiativeTransferTests/CollidingBeams_2D/run.sh b/examples/RadiativeTransferTests/CollidingBeams_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..e7af6a0ec80635dddb58397754590ca8854a732c
--- /dev/null
+++ b/examples/RadiativeTransferTests/CollidingBeams_2D/run.sh
@@ -0,0 +1,31 @@
+#!/bin/bash
+
+# exit if anything fails
+set -e
+set -o pipefail
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e glassPlane_128.hdf5 ]
+then
+ echo "Fetching initial glass file for the 2D RT advection example..."
+ ./getGlass.sh
+fi
+if [ ! -e collision_2D.hdf5 ]
+then
+ echo "Generating initial conditions for the 2D RT advection example..."
+ python3 makeIC.py
+fi
+
+# Run SWIFT with RT
+../../../swift \
+ --hydro \
+ --threads=4 \
+ --verbose=0 \
+ --radiation \
+ --stars \
+ --feedback \
+ --external-gravity \
+ ./rt_collision2D.yml 2>&1 | tee output.log
+
+python3 ./plotEnergies.py
+python3 ./plotSolution.py