diff --git a/Makefile.am b/Makefile.am
index 2c1af34d8df7b0bf83e6421c02db711c160cae3e..3ca9fd5e746dcbfa55d85e3632b3466d9ebc21ab 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -27,8 +27,10 @@ endif
SUBDIRS += src argparse examples doc tests tools
if HAVEEAGLECOOLING
SUBDIRS += examples/Cooling/CoolingRates
-endif
+DIST_SUBDIRS = $(SUBDIRS)
+else
DIST_SUBDIRS = $(SUBDIRS) examples/Cooling/CoolingRates
+endif
# Common flags
MYFLAGS =
diff --git a/configure.ac b/configure.ac
index 672dcd5680708de4abaf436be283e4e429e011a9..46271dc8699a41835cce4b04f10213a4e4056c48 100644
--- a/configure.ac
+++ b/configure.ac
@@ -2316,6 +2316,11 @@ case "$with_cooling" in
with_cooling_name=$with_cooling
;;
grackle_*)
+
+ if test "$have_grackle" != "yes"; then
+ AC_MSG_ERROR([Grackle cooling: You need the grackle library for Grackle cooling. (--with-grackle=PATH)])
+ fi
+
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
primordial_chemistry=${with_cooling#*_}
AC_DEFINE_UNQUOTED([COOLING_GRACKLE_MODE], [$primordial_chemistry], [Grackle chemistry network])
diff --git a/examples/Cooling/CoolingBox/coolingBox.yml b/examples/Cooling/CoolingBox/coolingBox.yml
index 9087a6087be4e4e5ea905b3fcfb85627bc02c413..7f4393781496ecc71e81de2544638dcda4d45c58 100644
--- a/examples/Cooling/CoolingBox/coolingBox.yml
+++ b/examples/Cooling/CoolingBox/coolingBox.yml
@@ -81,3 +81,13 @@ EAGLEChemistry: # Solar abundances
GEARPressureFloor:
jeans_factor: 0. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+
+COLIBRECooling:
+ dir_name: ./UV_dust1_CR1_G1_shield1.hdf5 # Location of the cooling tables
+ H_reion_z: 7.5 # Redshift of Hydrogen re-ionization (Planck 2018)
+ H_reion_eV_p_H: 2.0
+ He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
+ He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
+ He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
+ delta_logTEOS_subgrid_properties: 0.3 # delta log T above the EOS below which the subgrid properties use Teq assumption
+ rapid_cooling_threshold: 0.333333 # Switch to rapid cooling regime for dt / t_cool above this threshold.
diff --git a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py
index cb624be3cfb1f0f803cfce7a14fb1a772cccf515..2fe095ec948f78a8575561c62596bb078f3b67f0 100644
--- a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py
+++ b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py
@@ -36,6 +36,9 @@ def get_data_dump(metadata):
+ "$\\bf{Hydrodynamics}$\n"
+ metadata.hydro_info
+ "\n\n"
+ + "$\\bf{Cooling}$\n"
+ + metadata.subgrid_scheme["Cooling Model"].decode("utf-8")
+ + "\n\n"
+ "$\\bf{Viscosity}$\n"
+ viscosity
+ "\n\n"
@@ -68,8 +71,6 @@ def setup_axes():
for a in ax[:-1]:
a.set_xlim(0, 100)
- fig.tight_layout(pad=0.5)
-
return fig, ax
@@ -103,12 +104,15 @@ def get_data(handle: float, n_snaps: int):
try:
energies.append(
np.mean(
- (data.gas.internal_energies * data.gas.masses).to(erg).value
+ (data.gas.internal_energies * data.gas.masses)
+ .astype(np.float64)
+ .to(erg)
+ .value
)
* data.metadata.scale_factor ** (2)
)
radiated_energies.append(
- np.mean(data.gas.radiated_energy.to(erg).value)
+ np.mean(data.gas.radiated_energy.astype(np.float64).to(erg).value)
)
except AttributeError:
continue
@@ -152,7 +156,7 @@ def plot_single_data(
label_radiated = ""
label_sum = ""
- ax[2].plot(data[0], data[3], label=label_energy, ls="dotted", C=f"C{run}")
+ ax[2].plot(data[0], data[3], label=label_energy, ls="dotted", color=f"C{run}")
# ax[2].plot(data[0], data[4], label=label_radiated, ls="dashed", C=f"C{run}")
@@ -222,4 +226,5 @@ if __name__ == "__main__":
fig, ax = make_plot(handles, names, timestep_names)
+ fig.tight_layout(pad=0.5)
fig.savefig("redshift_dependence.png", dpi=300)
diff --git a/examples/Cooling/CoolingSedovBlast_3D/README b/examples/Cooling/CoolingSedovBlast_3D/README
new file mode 100644
index 0000000000000000000000000000000000000000..361cdd69eac8108a5603cc3c2501ae1ae425d912
--- /dev/null
+++ b/examples/Cooling/CoolingSedovBlast_3D/README
@@ -0,0 +1,29 @@
+Realistic SNe Sedov-Taylor explosion (realistic in the sense that it uses the
+right order of magnitude for all variables). Meant to be run with const Lambda
+cooling, but should also work with other cooling (the workflow also obtains
+the EAGLE cooling tables).
+
+The included Makefile (run.make) contains the full workflow that runs the
+simulation and creates a plot of the energy evolution, the evolution of the
+shock radius and velocity, and a movie of the density, velocity and temperature
+profile.
+
+To run the test, use
+ make -f run.make
+this will run the simulation in a subdirectory called "run".
+
+To run in a different subdirectory (useful for comparisons), use
+ make -f run.make folder=my_subdirectory_name
+
+The energy and shock evolution scripts (plot_energy.py and plot_time_profile.py)
+can also be run in comparison mode:
+ python3 plot_energy.py \
+ subdirectory1/statistics.txt subdirectory1/statistics.txt \
+ energy_1_vs_2.png \
+ --names "Run 1" "Run 2"
+and
+ python3 plot_time_evolution.py \
+ subdirectory1/profile.txt subdirectory1/profile.txt \
+ time_evolution_1_vs_2.png \
+ --names "Run 1" "Run 2"
+(there is not limit on the number of runs that can be included).
diff --git a/examples/Cooling/CoolingSedovBlast_3D/getGlass.sh b/examples/Cooling/CoolingSedovBlast_3D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..d5c5f590ac37c9c9431d626a2ea61b0c12c1513c
--- /dev/null
+++ b/examples/Cooling/CoolingSedovBlast_3D/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
diff --git a/examples/Cooling/CoolingSedovBlast_3D/get_time_profile.py b/examples/Cooling/CoolingSedovBlast_3D/get_time_profile.py
new file mode 100644
index 0000000000000000000000000000000000000000..00a5bb4817fc45088f4af60cdd95475d1d0e25cd
--- /dev/null
+++ b/examples/Cooling/CoolingSedovBlast_3D/get_time_profile.py
@@ -0,0 +1,69 @@
+import numpy as np
+import h5py
+import argparse
+import scipy.stats as stats
+import unyt
+
+argparser = argparse.ArgumentParser()
+argparser.add_argument("input", nargs="+")
+argparser.add_argument("output")
+args = argparser.parse_args()
+
+nfile = len(args.input)
+
+time = np.zeros(nfile)
+radius = np.zeros(nfile)
+velocity = np.zeros(nfile)
+
+for ifile, file in enumerate(args.input):
+ with h5py.File(file, "r") as handle:
+ coords = handle["PartType0/Coordinates"][:]
+ rho = handle["PartType0/Densities"][:]
+ vs = handle["PartType0/Velocities"][:]
+ t = handle["Header"].attrs["Time"][0]
+ box = handle["Header"].attrs["BoxSize"][:]
+
+ units = dict(handle["Units"].attrs)
+ uM = (units["Unit mass in cgs (U_M)"][0] * unyt.g).in_base("galactic")
+ uL = (units["Unit length in cgs (U_L)"][0] * unyt.cm).in_base("galactic")
+ ut = (units["Unit time in cgs (U_t)"][0] * unyt.s).in_base("galactic")
+
+ coords = (coords * uL).in_base("galactic")
+ rho = (rho * uM / uL ** 3).in_base("galactic")
+ vs = (vs * uL / ut).in_base("galactic")
+ t = (t * ut).in_base("galactic")
+ box = (box * uL).in_base("galactic")
+
+ coords -= 0.5 * box[None, :]
+
+ x = np.sqrt((coords ** 2).sum(axis=1))
+ v = (coords * vs).sum(axis=1) / x
+
+ x.convert_to_units("kpc")
+ v.convert_to_units("km/s")
+ t.convert_to_units("Myr")
+
+ rhohist, _, _ = stats.binned_statistic(x, rho, statistic="median", bins=100)
+ vhist, edges, _ = stats.binned_statistic(x, v, statistic="mean", bins=100)
+ mids = 0.5 * (edges[1:] + edges[:-1])
+
+ rhohist[np.isnan(rhohist)] = 0.0
+ imax = np.argmax(rhohist)
+
+ time[ifile] = t
+ radius[ifile] = mids[imax]
+ velocity[ifile] = vhist[imax]
+
+isort = np.argsort(time)
+
+time = time[isort]
+radius = radius[isort]
+velocity = velocity[isort]
+
+with open(args.output, "w") as handle:
+ handle.write("# Radial profile of shock wave\n")
+ handle.write("# Column 0: time (Myr)\n")
+ handle.write("# Column 1: radius (kpc)\n")
+ handle.write("# Column 2: velocity (km/s)\n")
+ for t, r, v in zip(time, radius, velocity):
+ handle.write(f"{t:.5e}\t{r:.5e}\t{v:.5e}\n")
diff --git a/examples/Cooling/CoolingSedovBlast_3D/makeIC.py b/examples/Cooling/CoolingSedovBlast_3D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..3449d73253422428a6c4367cc33d08f429560da9
--- /dev/null
+++ b/examples/Cooling/CoolingSedovBlast_3D/makeIC.py
@@ -0,0 +1,149 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+import h5py
+from numpy import *
+
+# Generates a swift IC file for the Sedov blast test in a periodic cubic box
+
+# Parameters
+gamma = 5.0 / 3.0 # Gas adiabatic index
+n0 = 0.1 # cm^-3
+mp = 1.67e-24 # g
+kB = 1.381e-16 # erg K^-1
+pc = 3.086e18 # cm
+rho0 = n0 * mp # Background density
+T0 = 1.0e4 # Background pressure
+E0 = 1.0e51 # Energy of the explosion
+N_inject = 32 # Number of particles in which to inject energy
+L = 256.0 * pc
+fileName = "sedov.hdf5"
+
+print(rho0, "cm s^-3")
+
+uL = 1.0e3 * pc
+uM = 1.99e30
+uv = 1.0e10
+ut = uL / uv
+uU = uv ** 2
+print("ut:", ut / 3.154e7, "yr")
+
+# ---------------------------------------------------
+glass = h5py.File("glassCube_64.hdf5", "r")
+
+# Read particle positions and h from the glass
+pos = glass["/PartType0/Coordinates"][:, :]
+h = glass["/PartType0/SmoothingLength"][:]
+
+pos *= L
+h *= L
+
+for i in range(3):
+ pos[pos[:, i] < 0.0, i] += L
+ pos[pos[:, i] >= L, i] -= L
+
+numPart = size(h)
+
+# newpos = zeros((numPart*8,3))
+# newh = zeros(numPart*8)
+# for ix in range(2):
+# for iy in range(2):
+# for iz in range(2):
+# offset = ix*4+iy*2+iz
+# newpos[offset*numPart:(offset+1)*numPart,:] = 0.5*pos[:,:]
+# newpos[offset*numPart:(offset+1)*numPart,0] += 0.5*ix*L
+# newpos[offset*numPart:(offset+1)*numPart,1] += 0.5*iy*L
+# newpos[offset*numPart:(offset+1)*numPart,2] += 0.5*iz*L
+# newh[offset*numPart:(offset+1)*numPart] = 0.5*h[:]
+
+# pos = newpos
+# h = newh
+print(h, h.min(), h.max())
+
+numPart = size(h)
+vol = L ** 3
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+r = zeros(numPart)
+
+r = sqrt(
+ (pos[:, 0] - 0.5 * L) ** 2 + (pos[:, 1] - 0.5 * L) ** 2 + (pos[:, 2] - 0.5 * L) ** 2
+)
+m[:] = rho0 * vol / numPart
+# u[:] = P0 / (rho0 * (gamma - 1))
+u[:] = kB * T0 / ((gamma - 1.0) * mp)
+
+print(u.mean(), E0 / (N_inject * m[0]))
+print(E0 / (N_inject * m[0]) / u.mean())
+
+# Make the central particle detonate
+index = argsort(r)
+u[index[0:N_inject]] = E0 / (N_inject * m[0])
+
+pos /= uL
+h /= uL
+L /= uL
+m /= uM
+u /= uU
+
+print("L:", L)
+print("m:", m.mean())
+print("h:", h.mean())
+print("u:", u.mean(), u.min(), u.max())
+print("pos:", pos.min(), pos.max())
+
+# --------------------------------------------------
+
+# File
+file = h5py.File(fileName, "w")
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [L, L, L]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 3
+
+# Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = uL
+grp.attrs["Unit mass in cgs (U_M)"] = uM
+grp.attrs["Unit time in cgs (U_t)"] = ut
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
+
+# Particle group
+grp = file.create_group("/PartType0")
+grp.create_dataset("Coordinates", data=pos, dtype="d")
+grp.create_dataset("Velocities", data=v, dtype="f")
+grp.create_dataset("Masses", data=m, dtype="f")
+grp.create_dataset("SmoothingLength", data=h, dtype="f")
+grp.create_dataset("InternalEnergy", data=u, dtype="f")
+grp.create_dataset("ParticleIDs", data=ids, dtype="L")
+
+file.close()
diff --git a/examples/Cooling/CoolingSedovBlast_3D/make_visualisations.make b/examples/Cooling/CoolingSedovBlast_3D/make_visualisations.make
new file mode 100644
index 0000000000000000000000000000000000000000..8c31a20fc3a0fcfbff7ce282bccd3a0dd452d1a3
--- /dev/null
+++ b/examples/Cooling/CoolingSedovBlast_3D/make_visualisations.make
@@ -0,0 +1,20 @@
+folder=.
+snaps=$(shell ls ${folder}/sedov_*.hdf5)
+imgs=$(patsubst ${folder}/sedov_%.hdf5,${folder}/SedovCooling_%.png,$(snaps))
+
+all: ${folder}/profile.png ${folder}/SedovCooling.mp4 ${folder}/energy.png
+
+${folder}/energy.png: ${folder}/statistics.txt
+ python3 plot_energy.py ${folder}/statistics.txt ${folder}/energy.png
+
+${folder}/profile.png: ${folder}/profile.txt
+ python3 plot_time_profile.py ${folder}/profile.txt ${folder}/profile.png
+
+${folder}/profile.txt: $(snaps)
+ python3 get_time_profile.py $(snaps) ${folder}/profile.txt
+
+${folder}/SedovCooling.mp4: $(imgs)
+ ffmpeg -y -framerate 10 -pattern_type glob -i "${folder}/SedovCooling_*.png" -c:v libx264 ${folder}/SedovCooling.mp4
+
+${folder}/SedovCooling_%.png: ${folder}/sedov_%.hdf5
+ python3 plot_profile.py $< $@
diff --git a/examples/Cooling/CoolingSedovBlast_3D/plot_energy.py b/examples/Cooling/CoolingSedovBlast_3D/plot_energy.py
new file mode 100644
index 0000000000000000000000000000000000000000..1c1e8a9d5394f3ae0270acbaa0210188f0100044
--- /dev/null
+++ b/examples/Cooling/CoolingSedovBlast_3D/plot_energy.py
@@ -0,0 +1,49 @@
+import numpy as np
+import argparse
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as pl
+
+argparser = argparse.ArgumentParser()
+argparser.add_argument("input", nargs="+")
+argparser.add_argument("output")
+argparser.add_argument("--names", "-n", nargs="+")
+args = argparser.parse_args()
+
+for ifile, file in enumerate(args.input):
+
+ data = np.loadtxt(
+ file,
+ usecols=(1, 13, 14, 16),
+ dtype=[
+ ("Time", np.float32),
+ ("Ekin", np.float32),
+ ("Etherm", np.float32),
+ ("Erad", np.float32),
+ ],
+ )
+
+ color = f"C{ifile}"
+
+ if args.names:
+ pl.plot([], [], "-", color=color, label=args.names[ifile])
+
+ pl.plot(data["Time"], data["Ekin"], "-.", color=color)
+ pl.plot(data["Time"], data["Etherm"], "--", color=color)
+ pl.plot(data["Time"], data["Erad"], ":", color=color)
+ pl.plot(
+ data["Time"], data["Ekin"] + data["Etherm"] + data["Erad"], "-", color=color
+ )
+
+pl.plot([], [], "k-", label="Total energy")
+pl.plot([], [], "k-.", label="Kinetic energy")
+pl.plot([], [], "k--", label="Thermal energy")
+pl.plot([], [], "k:", label="Radiated energy")
+pl.ylabel("Energy")
+pl.xlabel("Time")
+
+pl.legend(loc="best", ncol=3)
+
+pl.tight_layout()
+pl.savefig(args.output, dpi=300)
diff --git a/examples/Cooling/CoolingSedovBlast_3D/plot_profile.py b/examples/Cooling/CoolingSedovBlast_3D/plot_profile.py
new file mode 100644
index 0000000000000000000000000000000000000000..4c242da8bb3bc6cb2865cd792328c45686ad1db9
--- /dev/null
+++ b/examples/Cooling/CoolingSedovBlast_3D/plot_profile.py
@@ -0,0 +1,93 @@
+import numpy as np
+import h5py
+import argparse
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as pl
+import scipy.stats as stats
+import unyt
+
+
+def bin_vals(x, y, log=True):
+ stat = "median" if log else "mean"
+ hist, edges, _ = stats.binned_statistic(x, y, statistic=stat, bins=100)
+ mids = 0.5 * (edges[1:] + edges[:-1])
+ return mids, hist
+
+
+argparser = argparse.ArgumentParser()
+argparser.add_argument("input")
+argparser.add_argument("output")
+args = argparser.parse_args()
+
+x = None
+rho = None
+T = None
+v = None
+time = None
+with h5py.File(args.input, "r") as handle:
+ coords = handle["PartType0/Coordinates"][:]
+ rho = handle["PartType0/Densities"][:]
+ T = handle["PartType0/Temperatures"][:]
+ vs = handle["PartType0/Velocities"][:]
+ time = handle["Header"].attrs["Time"][0]
+ box = handle["Header"].attrs["BoxSize"][:]
+
+ units = dict(handle["Units"].attrs)
+ uM = (units["Unit mass in cgs (U_M)"][0] * unyt.g).in_base("galactic")
+ uL = (units["Unit length in cgs (U_L)"][0] * unyt.cm).in_base("galactic")
+ ut = (units["Unit time in cgs (U_t)"][0] * unyt.s).in_base("galactic")
+
+ coords = (coords * uL).in_base("galactic")
+ rho = (rho * uM / uL ** 3).in_base("galactic")
+ T = T * unyt.K
+ vs = (vs * uL / ut).in_base("galactic")
+ time = (time * ut).in_base("galactic")
+ box = (box * uL).in_base("galactic")
+
+ coords -= 0.5 * box[None, :]
+
+ x = np.sqrt((coords ** 2).sum(axis=1))
+ v = (coords * vs).sum(axis=1) / x
+
+rhohist, edges, _ = stats.binned_statistic(x, rho, statistic="median", bins=100)
+mids = 0.5 * (edges[1:] + edges[:-1])
+rhohist[np.isnan(rhohist)] = 0.0
+imax = np.argmax(rhohist)
+xmax = mids[imax]
+
+x.name = "radius"
+x.convert_to_units("kpc")
+v.name = "radial velocity"
+v.convert_to_units("km/s")
+rho.name = "density"
+rho.convert_to_units("g/cm**3")
+T.name = "temperature"
+
+fig, ax = pl.subplots(1, 3, figsize=(8, 4), sharex=True)
+
+with unyt.matplotlib_support:
+ ax[0].semilogy(x, rho, ".")
+ b, h = bin_vals(x, rho)
+ ax[0].semilogy(b, h, "--")
+
+ ax[1].plot(x, v, ".")
+ b, h = bin_vals(x, v, log=False)
+ ax[1].plot(b, h, "--")
+
+ ax[2].semilogy(x, T, ".")
+ b, h = bin_vals(x, T)
+ ax[2].semilogy(b, h, "--")
+
+for a in ax:
+ a.axvline(x=xmax, color="k", linestyle="--")
+
+ax[0].set_ylim(1.0e-28, 1.0e-24)
+ax[1].set_ylim(-10.0, 200.0)
+ax[2].set_ylim(10.0, 1.0e8)
+
+ax[1].set_title(f"t = {time:.2e}")
+
+pl.tight_layout()
+pl.savefig(args.output, dpi=300)
diff --git a/examples/Cooling/CoolingSedovBlast_3D/plot_time_profile.py b/examples/Cooling/CoolingSedovBlast_3D/plot_time_profile.py
new file mode 100644
index 0000000000000000000000000000000000000000..424850f3fef6544c9027733adb223d059cd3261e
--- /dev/null
+++ b/examples/Cooling/CoolingSedovBlast_3D/plot_time_profile.py
@@ -0,0 +1,41 @@
+import numpy as np
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as pl
+import argparse
+import unyt
+
+argparser = argparse.ArgumentParser()
+argparser.add_argument("input", nargs="+")
+argparser.add_argument("output")
+argparser.add_argument("--names", "-n", nargs="+")
+args = argparser.parse_args()
+
+fig, ax = pl.subplots(1, 2, sharex=True)
+
+for ifile, file in enumerate(args.input):
+ data = np.loadtxt(
+ file,
+ dtype=[("time", np.float32), ("radius", np.float32), ("velocity", np.float32)],
+ )
+
+ t = data["time"] * unyt.Myr
+ t.name = "time"
+ r = data["radius"] * unyt.kpc
+ r.name = "radius"
+ v = data["velocity"] * unyt.km / unyt.s
+ v.name = "velocity"
+
+ label = None
+ if args.names:
+ label = args.names[ifile]
+ with unyt.matplotlib_support:
+ ax[0].plot(t, r, "-", label=label)
+ ax[1].plot(t, v, "-", label=label)
+
+if args.names:
+ ax[1].legend(loc="best")
+
+pl.tight_layout()
+pl.savefig(args.output, dpi=300)
diff --git a/examples/Cooling/CoolingSedovBlast_3D/run.make b/examples/Cooling/CoolingSedovBlast_3D/run.make
new file mode 100644
index 0000000000000000000000000000000000000000..976962e80e071cdaea273d4ff86a2e986f3d101e
--- /dev/null
+++ b/examples/Cooling/CoolingSedovBlast_3D/run.make
@@ -0,0 +1,29 @@
+folder=run
+pwd=$(shell pwd)
+
+${folder}/profile.png: ${folder}/sedov_0100.hdf5
+ make -f make_visualisations.make -j 16 folder=${folder}
+
+${folder}/sedov_0100.hdf5: ${folder}/swift ${folder}/sedov.hdf5 ${folder}/coolingtables/redshifts.dat ${folder}/sedov.yml
+ cd ${folder}; ./swift --threads 8 --hydro --cooling --limiter sedov.yml
+
+${folder}/sedov.hdf5: sedov.hdf5
+ ln -s ${pwd}/sedov.hdf5 ${folder}/sedov.hdf5
+
+${folder}/coolingtables/redshifts.dat: coolingtables/redshifts.dat
+ ln -s ${pwd}/coolingtables ${folder}/coolingtables
+
+${folder}/swift: ../../../swift
+ mkdir -p ${folder}; ln -s ${pwd}/../../../swift ${folder}/swift
+
+${folder}/sedov.yml: sedov.yml
+ ln -s ${pwd}/sedov.yml ${folder}/sedov.yml
+
+sedov.hdf5: glassCube_64.hdf5
+ python3 makeIC.py
+
+glassCube_64.hdf5:
+ bash getGlass.sh
+
+coolingtables/redshifts.dat:
+ bash ../getEagleCoolingTable.sh
diff --git a/examples/Cooling/CoolingSedovBlast_3D/sedov.yml b/examples/Cooling/CoolingSedovBlast_3D/sedov.yml
new file mode 100644
index 0000000000000000000000000000000000000000..352686068fb9d5a34d3617ab869b1a7ecd92eebf
--- /dev/null
+++ b/examples/Cooling/CoolingSedovBlast_3D/sedov.yml
@@ -0,0 +1,99 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.99e33 # g
+ UnitLength_in_cgs: 3.086e21 # cm
+ UnitVelocity_in_cgs: 1.e5 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 1.e-3 # The end time of the simulation (in internal units).
+ dt_min: 1.e-13 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1.e-3 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sedov # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1.e-5 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-5 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 1.e2
+ max_ghost_iterations: 1000
+ h_max: 0.04 # Maximum smoothing length. Since the density in the post-shock region becomes very low, we need to make sure that this does not exceed 1/3 of the periodic box (including kernel gamma ~ 2)
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./sedov.hdf5 # The file to read
+ periodic: 1
+
+EAGLECooling:
+ H_reion_z: 11.5
+ H_reion_eV_p_H: 2.0
+ He_reion_z_centre: 3.5
+ He_reion_z_sigma: 0.5
+ He_reion_eV_p_H: 2.0
+ dir_name: ./coolingtables/
+
+EAGLEChemistry:
+ init_abundance_metal: 0.014 # Mass fraction in *all* metals
+ init_abundance_Hydrogen: 0.70649785 # Mass fraction in Hydrogen
+ init_abundance_Helium: 0.28055534 # Mass fraction in Helium
+ init_abundance_Carbon: 2.0665436e-3 # Mass fraction in Carbon
+ init_abundance_Nitrogen: 8.3562563e-4 # Mass fraction in Nitrogen
+ init_abundance_Oxygen: 5.4926244e-3 # Mass fraction in Oxygen
+ init_abundance_Neon: 1.4144605e-3 # Mass fraction in Neon
+ init_abundance_Magnesium: 5.907064e-4 # Mass fraction in Magnesium
+ init_abundance_Silicon: 6.825874e-4 # Mass fraction in Silicon
+ init_abundance_Iron: 1.1032152e-3 # Mass fraction in Iron
+
+EoS:
+ isothermal_internal_energy: 1240419161676.6465
+
+LambdaTTableCooling:
+ file_name: cooling.dat
+
+ConstCooling:
+ cooling_rate: 1.e6
+ min_energy: 0.01
+ cooling_tstep_mult: 2
+
+LambdaCooling:
+ lambda_nH2_cgs: 2.e-23
+ cooling_tstep_mult: 0.1
+ rapid_cooling: 1
+
+Scheduler:
+ max_top_level_cells: 12
+ cell_max_size: 8000000
+ cell_sub_size_pair_hydro: 256000000
+ cell_sub_size_self_hydro: 32000
+ cell_sub_size_pair_stars: 256000000
+ cell_sub_size_self_stars: 32000
+ cell_sub_size_pair_grav: 256000000
+ cell_sub_size_self_grav: 32000
+ cell_split_size: 400
+ cell_subdepth_diff_grav: 4
+ cell_extra_parts: 0
+ cell_extra_sparts: 100
+ cell_extra_gparts: 0
+ cell_extra_bparts: 0
+ cell_extra_sinks: 0
+ engine_max_parts_per_ghost: 1000
+ engine_max_sparts_per_ghost: 1000
+ engine_max_parts_per_cooling: 10000
+ nr_queues: 4
+ dependency_graph_frequency: 0
+ task_level_output_frequency: 0
+ tasks_per_cell: 100
+ links_per_tasks: 25
+ mpi_message_limit: 4
diff --git a/examples/EAGLE_DMO_low_z/EAGLE_DMO_100/eagle_100.yml b/examples/EAGLE_DMO_low_z/EAGLE_DMO_100/eagle_100.yml
index 9042d2d1c07b35adaa5b747d314c225a8a540b9e..dcecb9edd2476a87741c95341000424b6745e083 100644
--- a/examples/EAGLE_DMO_low_z/EAGLE_DMO_100/eagle_100.yml
+++ b/examples/EAGLE_DMO_low_z/EAGLE_DMO_100/eagle_100.yml
@@ -42,10 +42,10 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
- MAC: geometric
+ MAC: adaptive
+ epsilon_fmm: 0.001
theta_cr: 0.7 # Opening angle (Multipole acceptance criterion)
- use_tree_below_softening: 1
- mesh_side_length: 256
+ mesh_side_length: 1024
comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
diff --git a/examples/EAGLE_DMO_low_z/EAGLE_DMO_12/eagle_12.yml b/examples/EAGLE_DMO_low_z/EAGLE_DMO_12/eagle_12.yml
index bddb5e35baf4e1fed97b79254ba00d5811fac126..cd1165f6f4216f749d0164eededd036325b01bfb 100644
--- a/examples/EAGLE_DMO_low_z/EAGLE_DMO_12/eagle_12.yml
+++ b/examples/EAGLE_DMO_low_z/EAGLE_DMO_12/eagle_12.yml
@@ -42,10 +42,10 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
- MAC: geometric
+ MAC: adaptive
+ epsilon_fmm: 0.001
theta_cr: 0.7 # Opening angle (Multipole acceptance criterion)
- use_tree_below_softening: 1
- mesh_side_length: 32
+ mesh_side_length: 128
comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
diff --git a/examples/EAGLE_DMO_low_z/EAGLE_DMO_25/eagle_25.yml b/examples/EAGLE_DMO_low_z/EAGLE_DMO_25/eagle_25.yml
index 4fd5d226735ba6c49055de9b381b908c0beb5c52..163aba63e7cd7c03a8cc19f55eda3015fdf0a322 100644
--- a/examples/EAGLE_DMO_low_z/EAGLE_DMO_25/eagle_25.yml
+++ b/examples/EAGLE_DMO_low_z/EAGLE_DMO_25/eagle_25.yml
@@ -42,10 +42,10 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
- MAC: geometric
+ MAC: adaptive
+ epsilon_fmm: 0.001
theta_cr: 0.7 # Opening angle (Multipole acceptance criterion)
- use_tree_below_softening: 1
- mesh_side_length: 64
+ mesh_side_length: 256
comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
diff --git a/examples/EAGLE_DMO_low_z/EAGLE_DMO_50/eagle_50.yml b/examples/EAGLE_DMO_low_z/EAGLE_DMO_50/eagle_50.yml
index 6f19042e6850b6960f03817116f9c2ea73c70bf5..b77a2bcabd43d94c21846fd40be7b1d499b292bb 100644
--- a/examples/EAGLE_DMO_low_z/EAGLE_DMO_50/eagle_50.yml
+++ b/examples/EAGLE_DMO_low_z/EAGLE_DMO_50/eagle_50.yml
@@ -41,10 +41,10 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
- MAC: geometric
+ MAC: adaptive
+ epsilon_fmm: 0.001
theta_cr: 0.7 # Opening angle (Multipole acceptance criterion)
- use_tree_below_softening: 1
- mesh_side_length: 128
+ mesh_side_length: 512
comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
diff --git a/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/README b/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/README
new file mode 100644
index 0000000000000000000000000000000000000000..3cf3226b923b9b9206377770fac9fb1f6f68f06c
--- /dev/null
+++ b/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/README
@@ -0,0 +1,45 @@
+1D advection test for radiative transfer.
+This test is identical to the examples/RadiativeTransferTests/Advection_1D
+with the difference that particles aren't equally spaced on a line, but
+separated into three regions with decreasing intervals between particles
+in each of the region.
+
+More concretely, for a box of size 1, the separation is given as
+
+dx if x < 0.33
+dx/2 if 0.33 < x < 0.66
+dx/4 if x >= 0.66
+
+This leads to the three regions having different time step sizes for RT.
+
+
+Test that your method is TVD and the propagation speed of the photons is
+correct. The ICs set up three photon groups:
+- The first is a top hat function initial distribution where outside values
+ are zero.
+- The second is a top hat function initial distribution where outside values
+ are nonzero. This distinction is important to test because photon energies
+ can't be negative, so these cases need to be tested individually.
+- the third is a smooth Gaussian.
+
+This way, you can test multiple initial condition scenarios simultaneously.
+There are no stars to act as sources. Also make sure that you choose your
+photon frequencies in a way that doesn't interact with gas!
+
+The ICs are created to be compatible with GEAR_RT and SPHM1RT. Recommended configuration:
+GEAR_RT:
+ --with-rt=GEAR_3 --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
+
+SPHM1RT:
+ --with-rt=SPHM1RT_3 --with-hydro-dimension=1 --with-stars=basic --with-sundials=$SUNDIALS_ROOT
+
+IMPORTANT: Need SUNDIALS version = 5 .
+SUNDIALS_ROOT is the root directory that contains the lib and include directories, e.g. on cosma:
+SUNDIALS_ROOT=/cosma/local/sundials/5.1.0/
+
+Note that in SPHM1RT any (SPH) hydro scheme is compatible.
+
+Note that if you want to use a reduced speed of light for this test, you also
+need to adapt the fluxes in the initial conditions! They are generated assuming
+that the speed of light is not reduced.
diff --git a/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/makeIC.py b/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/makeIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..8d61caffb0f9df90fa55c6addab1a4c078162b41
--- /dev/null
+++ b/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/makeIC.py
@@ -0,0 +1,207 @@
+#!/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
+dx_init = boxsize / 1000
+
+# filename of ICs to be generated
+outputfilename = "advection_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:
+ # -------------------
+
+ if x[0] < 0.33 * boxsize:
+ E = 0.0
+ elif x[0] < 0.66 * boxsize:
+ E = 1.0
+ else:
+ E = 0.0
+
+ # Assuming all photons flow in only one direction
+ # (optically thin regime, "free streaming limit"),
+ # we have that |F| = c * E
+ F = np.zeros(3, dtype=np.float64)
+ F[0] = unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
+
+ E_list.append(E)
+ F_list.append(F)
+
+ # Group 2 Photons:
+ # -------------------
+
+ if x[0] < 0.33 * boxsize:
+ E = 1.0
+ elif x[0] < 0.66 * boxsize:
+ E = 3.0
+ else:
+ E = 1.0
+
+ F = np.zeros(3, dtype=np.float64)
+ F[0] = unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
+
+ E_list.append(E)
+ F_list.append(F)
+
+ # Group 3 Photons:
+ # -------------------
+ sigma = 0.1 * boxsize
+ mean = 0.5 * boxsize
+ amplitude = 2.0
+
+ E = amplitude * np.exp(-(x[0] - mean) ** 2 / (2 * sigma ** 2))
+ F = np.zeros(3, dtype=np.float64)
+ F[0] = 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_list = []
+ volumes_list = []
+
+ x = 0.5 * dx_init
+ while x < boxsize:
+ if x < 0.33 * boxsize:
+ dx = dx_init
+ elif 0.33 * boxsize <= x <= 0.66 * boxsize:
+ dx = dx_init / 2
+ else:
+ dx = dx_init / 4
+ x = x + dx
+ xp_list.append(x)
+ volumes_list.append(dx)
+
+ n_p = len(xp_list)
+ xp = unyt.unyt_array(np.zeros((n_p, 3), dtype=np.float64), boxsize.units)
+ volumes = unyt.unyt_array(np.zeros((n_p), dtype=np.float64), boxsize.units ** 3)
+ for i in range(n_p):
+ xp[i, 0] = xp_list[i]
+ volumes[i] = volumes_list[i]
+
+ 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] / volumes[p]
+ # Fsetname = "PhotonFluxesGroup{0:d}".format(g + 1)
+ # parts[Fsetname][p] = Flux[g] / volumes[p]
+ # Esetname = "PhotonEnergiesGroup{0:d}".format(g + 1)
+ # parts[Esetname][p] = E[g]
+ # Fsetname = "PhotonFluxesGroup{0:d}".format(g + 1)
+ # parts[Fsetname][p] = Flux[g]
+ Esetname = "PhotonEnergiesGroup{0:d}".format(g + 1)
+ parts[Esetname][p] = E[g] * volumes[p]
+ Fsetname = "PhotonFluxesGroup{0:d}".format(g + 1)
+ parts[Fsetname][p] = Flux[g] * volumes[p]
+
+ # 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/AdvectionDifferentTimeStepSizes_1D/plotSolution.py b/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/plotSolution.py
new file mode 100755
index 0000000000000000000000000000000000000000..fb6a962976e3e94f1fa37fb9dbe88d889f6eabd8
--- /dev/null
+++ b/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/plotSolution.py
@@ -0,0 +1,231 @@
+#!/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 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
+fancy = True # fancy up the plots a bit
+plot_analytical_solutions = True # overplot analytical solution
+
+# properties for all scatterplots
+scatterplot_kwargs = {
+ "facecolor": "red",
+ "s": 4,
+ "alpha": 0.6,
+ "linewidth": 0.0,
+ "marker": ".",
+}
+
+# properties for all analytical solution plots
+analytical_solution_kwargs = {"linewidth": 1.0, "ls": "--", "c": "k", "alpha": 0.5}
+
+# -----------------------------------------------------------------------
+
+if plot_analytical_solutions:
+ from makeIC import initial_condition
+
+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"))
+
+ boxsize = meta.boxsize[0]
+ 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()
+
+ # get analytical solutions
+ if plot_analytical_solutions:
+
+ time = meta.time
+ speed = meta.reduced_lightspeed
+
+ advected_positions = data.gas.coordinates[:].copy()
+ advected_positions[:, 0] -= speed * time
+ nparts = advected_positions.shape[0]
+ # add periodicity corrections
+ negatives = advected_positions < 0.0
+ if negatives.any():
+ while advected_positions.min() < 0.0:
+ advected_positions[negatives] += boxsize
+ overshooters = advected_positions > boxsize
+ if overshooters.any():
+ while advected_positions.max() > boxsize:
+ advected_positions[overshooters] -= boxsize
+
+ analytical_solutions = np.zeros((nparts, ngroups), dtype=np.float64)
+ for p in range(part_positions.shape[0]):
+ E, F = initial_condition(advected_positions[p])
+ for g in range(ngroups):
+ analytical_solutions[p, g] = E[g]
+
+ fig = plt.figure(figsize=(5.05 * ngroups, 5.4), dpi=200)
+ figname = filename[:-5] + ".png"
+
+ for g in range(ngroups):
+
+ # plot energy density
+ new_attribute_str = "radiation_energy" + str(g + 1)
+ photon_energy = getattr(data.gas, new_attribute_str)
+
+ volumes = data.gas.masses / data.gas.densities
+ photon_energy_density = photon_energy / volumes
+
+ ax = fig.add_subplot(2, ngroups, g + 1)
+ s = np.argsort(part_positions)
+ if plot_analytical_solutions:
+ ax.plot(
+ part_positions[s],
+ analytical_solutions[s, g],
+ **analytical_solution_kwargs,
+ label="analytical solution",
+ )
+ ax.scatter(
+ part_positions,
+ photon_energy_density,
+ **scatterplot_kwargs,
+ label="simulation",
+ )
+ ax.legend()
+
+ ax.set_title("Group {0:2d}".format(g + 1))
+ if g == 0:
+ ax.set_ylabel(
+ "Energy Density [$"
+ + photon_energy_density.units.latex_representation()
+ + "$]"
+ )
+ ax.set_xlabel("x [$" + part_positions.units.latex_representation() + "$]")
+
+ # plot flux X
+ new_attribute_str = "radiation_flux" + str(g + 1) + "X"
+ photon_flux = getattr(data.gas, new_attribute_str)
+
+ if scheme.startswith("SPH M1closure"):
+ photon_flux = photon_flux / volumes
+
+ photon_flux = photon_flux.to("erg/cm**2/s")
+
+ ax2 = fig.add_subplot(2, ngroups, g + 1 + ngroups)
+ ax2.scatter(part_positions, photon_flux, **scatterplot_kwargs)
+
+ if g == 0:
+ ax2.set_ylabel(
+ "Flux X [$" + photon_flux.units.latex_representation() + "$]"
+ )
+ ax2.set_xlabel("x [$" + part_positions.units.latex_representation() + "$]")
+
+ # 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
+
+
+if __name__ == "__main__":
+
+ snaplist = get_snapshot_list(snapshot_base)
+
+ for f in snaplist:
+ plot_photons(f)
diff --git a/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/rt_advection1D.yml b/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/rt_advection1D.yml
new file mode 100644
index 0000000000000000000000000000000000000000..d34aabc03d321e595a3e882c852ebaf7e4539deb
--- /dev/null
+++ b/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/rt_advection1D.yml
@@ -0,0 +1,81 @@
+MetaData:
+ run_name: RT_advection-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: 4.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: ./advection_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.4 # 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
+
+SPHM1RT:
+ cred: 2.99792458e10 # reduce the speed of light in the code unit
+ CFL_condition: 0.1 # CFL condition for RT, independent of hydro
+ photon_groups_Hz: [1., 2.] # Photon frequency group bin edges in Hz. Needs to be 1 less than the number of groups (N) requested during the configuration (--with-RT=SPHM1RT_N).
+ use_const_emission_rates: 1 # (Optional) use constant emission rates for stars as defined with star_emission_rates parameter
+ star_emission_rates: [1e-32, 1e-32, 1e-32] # (Optional) constant star emission rates for each photon frequency group to use if use_constant_emission_rates is set.
+ 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 # (Optional) skip the thermochemistry. This is intended only for debugging and testing the radiation transport, as it breaks the purpose of RT.
+ init_mass_fraction_metal: 0. # (Optional) Inital mass fraction of particle mass in *all* metals (if it is set, the initial fraction will be over-written.)
+ init_mass_fraction_Hydrogen: 0.752 # (Conditional) (if init_mass_fraction_metal != -1.0f) Inital mass fraction of particle mass in Hydrogen
+ init_mass_fraction_Helium: 0.248 # (Conditional) (if init_mass_fraction_metal != -1.0f) Inital mass fraction of particle mass in Helium
+ useabundances: 0 # (Optional) use the species abundances below, instead of reading from initial condition
+ init_species_abundance_e: 0.0 # (Conditional) (if useabundances==1) free electron abundances (in unit hydrogen number density:nH)
+ init_species_abundance_HI: 1.0 # (Conditional) (if useabundances==1) HI abundances (in unit hydrogen number density:nH)
+ init_species_abundance_HII: 0.0 # (Conditional) (if useabundances==1) HII abundances (in unit hydrogen number density:nH)
+ init_species_abundance_HeI: 0.08244680851 # (Conditional) (if useabundances==1) HeI abundances (in unit hydrogen number density:nH)
+ init_species_abundance_HeII: 0.0 # (Conditional) (if useabundances==1) HeII abundances (in unit hydrogen number density:nH)
+ init_species_abundance_HeIII: 0.0 # (Conditional) (if useabundances==1) HeIII abundances (in unit hydrogen number density:nH)
diff --git a/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/run.sh b/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..13e8331da60a1f432ff5c966c4f2da997bfff19a
--- /dev/null
+++ b/examples/RadiativeTransferTests/AdvectionDifferentTimeStepSizes_1D/run.sh
@@ -0,0 +1,23 @@
+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -f advection_1D.hdf5 ]; then
+ echo "Generating ICs"
+ python3 makeIC.py
+fi
+
+# Run SWIFT with RT
+../../../swift \
+ --hydro \
+ --threads=4 \
+ --verbose=0 \
+ --radiation \
+ --stars \
+ --feedback \
+ --external-gravity \
+ ./rt_advection1D.yml 2>&1 | tee output.log
+
+python3 ./plotSolution.py
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
diff --git a/src/cell.h b/src/cell.h
index c992171d714deded60e97e7df6205973f25cb0c4..5c407990a4de577a27d14908c9500e85b7c2251a 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -480,7 +480,10 @@ struct cell {
float dmin;
/*! ID of the previous owner, e.g. runner. */
- int owner;
+ short int owner;
+
+ /*! ID of a threadpool thread that maybe associated with this cell. */
+ short int tpid;
/*! ID of the node this cell lives on. */
int nodeID;
diff --git a/src/cell_unskip.c b/src/cell_unskip.c
index 538eeec0c5f01c9788447e713f709e4890bbf913..93d98b1b1910db988ed6372c66df711a32f9addd 100644
--- a/src/cell_unskip.c
+++ b/src/cell_unskip.c
@@ -3152,6 +3152,18 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
if (ci_active) {
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_transport);
}
+ } else if (ci_active) {
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+ /* If the local cell is inactive and the remote cell is active, we
+ * still need to receive stuff to be able to do the force interaction
+ * on this node as well.
+ * The gradient recv is only necessary in normal steps in case we need
+ * to sort, not during sub-cycles. */
+ if (!sub_cycle)
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_gradient);
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_transport);
+ if (sub_cycle) cell_set_skip_rt_sort_flag_up(ci);
+#endif
}
/* Is the foreign cell active and will need stuff from us? */
@@ -3164,6 +3176,19 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
scheduler_activate_send(s, cj->mpi.send, task_subtype_rt_transport,
ci_nodeID);
}
+ } else if (cj_active) {
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+ /* If the foreign cell is inactive, but the local cell is active,
+ * we still need to send stuff to be able to do the force interaction
+ * on both nodes.
+ * The gradient send is only necessary in normal steps in case we need
+ * to sort, not during sub-cycles. */
+ if (!sub_cycle)
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_rt_gradient,
+ ci_nodeID);
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_rt_transport,
+ ci_nodeID);
+#endif
}
} else if (cj_nodeID != nodeID) {
@@ -3180,6 +3205,18 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
if (cj_active) {
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_transport);
}
+ } else if (cj_active) {
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+ /* If the local cell is inactive and the remote cell is active, we
+ * still need to receive stuff to be able to do the force interaction
+ * on this node as well.
+ * The gradient recv is only necessary in normal steps in case we need
+ * to sort, not during sub-cycles. */
+ if (!sub_cycle)
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_gradient);
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_transport);
+ if (sub_cycle) cell_set_skip_rt_sort_flag_up(cj);
+#endif
}
/* Is the foreign cell active and will need stuff from us? */
@@ -3192,6 +3229,19 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
scheduler_activate_send(s, ci->mpi.send, task_subtype_rt_transport,
cj_nodeID);
}
+ } else if (ci_active) {
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+ /* If the foreign cell is inactive, but the local cell is active,
+ * we still need to send stuff to be able to do the force interaction
+ * on both nodes
+ * The gradient send is only necessary in normal steps in case we need
+ * to sort, not during sub-cycles. */
+ if (!sub_cycle)
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_rt_gradient,
+ cj_nodeID);
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_rt_transport,
+ cj_nodeID);
+#endif
}
}
#endif
@@ -3240,7 +3290,52 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
scheduler_activate(s, c->rt.rt_transport_out);
if (c->rt.rt_tchem != NULL) scheduler_activate(s, c->rt.rt_tchem);
if (c->rt.rt_out != NULL) scheduler_activate(s, c->rt.rt_out);
+ } else {
+#if defined(MPI_SYMMETRIC_FORCE_INTERACTION) && defined(WITH_MPI)
+ /* Additionally unskip force interactions between inactive local cell and
+ * active remote cell. (The cell unskip will only be called for active
+ * cells, so, we have to do this now, from the active remote cell). */
+ for (struct link *l = c->rt.rt_transport; l != NULL; l = l->next) {
+ struct task *t = l->t;
+ if (t->type != task_type_pair && t->type != task_type_sub_pair)
+ continue;
+
+ struct cell *ci = l->t->ci;
+ struct cell *cj = l->t->cj;
+
+ const int ci_active = cell_is_rt_active(ci, e);
+ const int cj_active = cell_is_rt_active(cj, e);
+ const int ci_nodeID = ci->nodeID;
+ const int cj_nodeID = cj->nodeID;
+ if ((!ci_active && ci_nodeID == nodeID && cj_active &&
+ cj_nodeID != nodeID) ||
+ (!cj_active && cj_nodeID == nodeID && ci_active &&
+ ci_nodeID != nodeID)) {
+ scheduler_activate(s, l->t);
+
+ if (!sub_cycle) {
+ /* Activate sorts only during main/normal steps. */
+ if (t->type == task_type_pair) {
+ atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
+ atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
+ ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
+ cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
+
+ /* Check the sorts and activate them if needed. */
+ cell_activate_rt_sorts(ci, t->flags, s);
+ cell_activate_rt_sorts(cj, t->flags, s);
+ }
+
+ /* Store current values of dx_max and h_max. */
+ else if (t->type == task_type_sub_pair) {
+ cell_activate_subcell_rt_tasks(ci, cj, s, sub_cycle);
+ }
+ }
+ }
+ }
+#endif
}
+
/* The rt_advance_cell_time tasks also run on foreign cells */
if (c->super != NULL && c->super->rt.rt_advance_cell_time != NULL)
scheduler_activate(s, c->super->rt.rt_advance_cell_time);
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index b09e6610d58b5ae531b232b7539afb32082d658a..64a2d35040e25032961d906c38d02e0f9325ef2f 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -178,12 +178,24 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
total_du_dt = -u_old_com / ((2.5f + 0.0001f) * dt_therm);
}
- /* Update the internal energy time derivative */
- hydro_set_comoving_internal_energy_dt(p, total_du_dt);
+ if (cooling->rapid_cooling) {
+ const float u_new_com = u_old_com + total_du_dt * dt_therm;
+ const float u_new_phys = u_new_com * cosmo->a_factor_internal_energy;
+ hydro_set_physical_internal_energy(p, xp, cosmo, u_new_phys);
+ hydro_set_drifted_physical_internal_energy(p, cosmo, u_new_phys);
+ hydro_set_physical_internal_energy_dt(p, cosmo, 0.);
+ } else {
+ /* Update the internal energy time derivative */
+ hydro_set_comoving_internal_energy_dt(p, total_du_dt);
+ }
+ const float actual_cooling_du_dt = total_du_dt - hydro_du_dt_com;
+ const float actual_cooling_du_dt_physical = actual_cooling_du_dt / cosmo->a /
+ cosmo->a *
+ cosmo->a_factor_internal_energy;
/* Store the radiated energy (assuming dt will not change) */
xp->cooling_data.radiated_energy +=
- -hydro_get_mass(p) * cooling_du_dt_physical * dt;
+ -hydro_get_mass(p) * actual_cooling_du_dt_physical * dt;
}
/**
@@ -415,6 +427,8 @@ static INLINE void cooling_init_backend(struct swift_params* parameter_file,
parser_get_param_double(parameter_file, "LambdaCooling:lambda_nH2_cgs");
cooling->cooling_tstep_mult = parser_get_opt_param_float(
parameter_file, "LambdaCooling:cooling_tstep_mult", FLT_MAX);
+ cooling->rapid_cooling = parser_get_opt_param_int(
+ parameter_file, "LambdaCooling:rapid_cooling", 0);
/* Some useful conversion values */
cooling->conv_factor_density_to_cgs =
@@ -463,6 +477,12 @@ static INLINE void cooling_print_backend(
else
message("Cooling function time-step size limited to %f of u/(du/dt)",
cooling->cooling_tstep_mult);
+
+ if (cooling->rapid_cooling) {
+ message("Using rapid cooling");
+ } else {
+ message("Using normal cooling");
+ }
}
/**
diff --git a/src/cooling/const_lambda/cooling_io.h b/src/cooling/const_lambda/cooling_io.h
index 0efba80c294b732e2e4e59d824143b879757e37d..9494731d047b4f1b7b03d03e85d128c286351ae9 100644
--- a/src/cooling/const_lambda/cooling_io.h
+++ b/src/cooling/const_lambda/cooling_io.h
@@ -51,6 +51,7 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour(
io_write_attribute_s(h_grp, "Cooling Model", "Constant Lambda");
io_write_attribute_d(h_grp, "Lambda/n_H^2 [cgs]", cooling->lambda_nH2_cgs);
+ io_write_attribute_i(h_grp, "Rapid cooling", cooling->rapid_cooling);
}
#endif
@@ -76,7 +77,7 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
struct io_props* list) {
list[0] = io_make_output_field_convert_part(
- "Temperature", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts, xparts,
+ "Temperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts, xparts,
convert_part_T, "Temperatures of the gas particles");
list[1] = io_make_output_field(
diff --git a/src/cooling/const_lambda/cooling_properties.h b/src/cooling/const_lambda/cooling_properties.h
index 8417c9112645713a43c1edf0600db98f604d030e..ecb5c9c8e25d41a1ada04727ad403d6db5bf3ce4 100644
--- a/src/cooling/const_lambda/cooling_properties.h
+++ b/src/cooling/const_lambda/cooling_properties.h
@@ -48,6 +48,9 @@ struct cooling_function_data {
/*! Constant multiplication factor for time-step criterion */
float cooling_tstep_mult;
+
+ /*! Use rapid cooling? */
+ int rapid_cooling;
};
#endif /* SWIFT_COOLING_PROPERTIES_CONST_LAMBDA_H */
diff --git a/src/cosmology.c b/src/cosmology.c
index d7dc6a3f9fae48fb72d17d098e88b4b2057eac8d..6bb6d4c8f9feeb281104b6ff7f0dcab2dbfb0a53 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -1509,7 +1509,7 @@ void cosmology_write_model(hid_t h_grp, const struct cosmology *c) {
io_write_attribute_d(h_grp, "T_nu_0 [eV]", c->T_nu_0_eV);
io_write_attribute_d(h_grp, "N_eff", c->N_eff);
io_write_attribute_d(h_grp, "N_ur", c->N_ur);
- io_write_attribute_d(h_grp, "N_nu", c->N_nu);
+ io_write_attribute_i(h_grp, "N_nu", c->N_nu);
if (c->N_nu > 0) {
io_write_attribute(h_grp, "M_nu_eV", DOUBLE, c->M_nu_eV, c->N_nu);
io_write_attribute(h_grp, "deg_nu", DOUBLE, c->deg_nu, c->N_nu);
diff --git a/src/engine_config.c b/src/engine_config.c
index b444c9a8ed357eb3e88728f6e2ff18e8bc049e20..3834398f812005087277640aabb8e493c47a5d4a 100644
--- a/src/engine_config.c
+++ b/src/engine_config.c
@@ -23,6 +23,7 @@
#include <config.h>
/* System includes. */
+#include <stdlib.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>
@@ -40,6 +41,7 @@
#include "part.h"
#include "pressure_floor.h"
#include "proxy.h"
+#include "rt.h"
#include "star_formation.h"
#include "star_formation_logger.h"
#include "stars_io.h"
@@ -195,6 +197,9 @@ void engine_config(int restart, int fof, struct engine *e,
e->restart_dt = 0;
e->run_fof = 0;
+ /* Seed rand(). */
+ srand(clocks_random_seed());
+
/* Allow repartitioning to be changed between restarts. On restart this is
* already allocated and freed on exit, so we need to copy over. */
#ifdef WITH_MPI
@@ -807,6 +812,12 @@ void engine_config(int restart, int fof, struct engine *e,
if (e->nodeID == 0)
message("Using %d threads in the thread-pool", nr_pool_threads);
+ /* Cells per thread buffer. */
+ e->s->cells_sub =
+ (struct cell **)calloc(nr_pool_threads + 1, sizeof(struct cell *));
+ e->s->multipoles_sub = (struct gravity_tensors **)calloc(
+ nr_pool_threads + 1, sizeof(struct gravity_tensors *));
+
/* First of all, init the barrier and lock it. */
if (swift_barrier_init(&e->wait_barrier, NULL, e->nr_threads + 1) != 0 ||
swift_barrier_init(&e->run_barrier, NULL, e->nr_threads + 1) != 0)
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index 7a2f38e4aa32abc062d39a906ed2d4bf4f69407c..edbcee18fd4015595de2cc8d017184d0f494b598 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -2813,6 +2813,32 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
sched, task_type_pair, task_subtype_rt_gradient, flags, 0, ci, cj);
t_rt_transport = scheduler_addtask(
sched, task_type_pair, task_subtype_rt_transport, flags, 0, ci, cj);
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+ /* The order of operations for an inactive local cell interacting
+ * with an active foreign cell is not guaranteed because the gradient
+ * iact loops don't exist in that case. So we need an explicit
+ * dependency here to have sorted cells. */
+
+ /* Make all force tasks depend on the sorts */
+ if (ci->hydro.super->rt.rt_sorts != NULL)
+ scheduler_addunlock(sched, ci->hydro.super->rt.rt_sorts,
+ t_rt_transport);
+ if (ci->hydro.super != cj->hydro.super) {
+ if (cj->hydro.super->rt.rt_sorts != NULL)
+ scheduler_addunlock(sched, cj->hydro.super->rt.rt_sorts,
+ t_rt_transport);
+ }
+ /* We need to ensure that a local inactive cell is sorted before
+ * the interaction in the transport loop. Local cells don't have an
+ * rt_sorts task. */
+ if (ci->hydro.super->hydro.sorts != NULL)
+ scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
+ t_rt_transport);
+ if ((ci->hydro.super != cj->hydro.super) &&
+ (cj->hydro.super->hydro.sorts != NULL))
+ scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
+ t_rt_transport);
+#endif
}
engine_addlink(e, &ci->hydro.force, t_force);
@@ -3615,6 +3641,32 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
t_rt_transport =
scheduler_addtask(sched, task_type_sub_pair,
task_subtype_rt_transport, flags, 0, ci, cj);
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+ /* The order of operations for an inactive local cell interacting
+ * with an active foreign cell is not guaranteed because the gradient
+ * iact loops don't exist in that case. So we need an explicit
+ * dependency here to have sorted cells. */
+
+ /* Make all force tasks depend on the sorts */
+ if (ci->hydro.super->rt.rt_sorts != NULL)
+ scheduler_addunlock(sched, ci->hydro.super->rt.rt_sorts,
+ t_rt_transport);
+ if (ci->hydro.super != cj->hydro.super) {
+ if (cj->hydro.super->rt.rt_sorts != NULL)
+ scheduler_addunlock(sched, cj->hydro.super->rt.rt_sorts,
+ t_rt_transport);
+ }
+ /* We need to ensure that a local inactive cell is sorted before
+ * the interaction in the transport loop. Local cells don't have
+ * an rt_sort task. */
+ if (ci->hydro.super->hydro.sorts != NULL)
+ scheduler_addunlock(sched, ci->hydro.super->hydro.sorts,
+ t_rt_transport);
+ if ((ci->hydro.super != cj->hydro.super) &&
+ (cj->hydro.super->hydro.sorts != NULL))
+ scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
+ t_rt_transport);
+#endif
}
engine_addlink(e, &ci->hydro.force, t_force);
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index e7f1d0263d3b13200a56886583f60c8cf196910a..b027156ff0371abccdfdf9d19e0b4cf1942e3076 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -903,6 +903,37 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
with_timestep_limiter);
}
}
+ }
+
+ /* Pair tasks between inactive local cells and active remote cells. */
+ if ((ci_nodeID != nodeID && cj_nodeID == nodeID && ci_active_rt &&
+ !cj_active_rt) ||
+ (ci_nodeID == nodeID && cj_nodeID != nodeID && !ci_active_rt &&
+ cj_active_rt)) {
+
+ if (t_subtype == task_subtype_rt_transport) {
+
+ scheduler_activate(s, t);
+
+ /* Set the correct sorting flags */
+ if (t_type == task_type_pair) {
+
+ /* Store some values. */
+ atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
+ atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
+ ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
+ cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
+
+ /* Check the sorts and activate them if needed. */
+ cell_activate_rt_sorts(ci, t->flags, s);
+ cell_activate_rt_sorts(cj, t->flags, s);
+ }
+
+ /* Store current values of dx_max and h_max. */
+ else if (t_type == task_type_sub_pair) {
+ cell_activate_subcell_rt_tasks(ci, cj, s, /*sub_cycle=*/0);
+ }
+ }
#endif
}
@@ -1389,6 +1420,14 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
scheduler_activate_recv(s, ci->mpi.recv,
task_subtype_rt_transport);
}
+ } else if (ci_active_rt) {
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+ /* If the local cell is inactive and the remote cell is active, we
+ * still need to receive stuff to be able to do the force
+ * interaction on this node as well. */
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_gradient);
+ scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_transport);
+#endif
}
/* Is the foreign cell active and will need stuff from us? */
@@ -1401,6 +1440,16 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
scheduler_activate_send(s, cj->mpi.send,
task_subtype_rt_transport, ci_nodeID);
}
+ } else if (cj_active_rt) {
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+ /* If the foreign cell is inactive, but the local cell is active,
+ * we still need to send stuff to be able to do the force
+ * interaction on both nodes */
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_rt_gradient,
+ ci_nodeID);
+ scheduler_activate_send(s, cj->mpi.send, task_subtype_rt_transport,
+ ci_nodeID);
+#endif
}
} else if (cj_nodeID != nodeID) {
@@ -1416,6 +1465,14 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
scheduler_activate_recv(s, cj->mpi.recv,
task_subtype_rt_transport);
}
+ } else if (cj_active_rt) {
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+ /* If the local cell is inactive and the remote cell is active, we
+ * still need to receive stuff to be able to do the force
+ * interaction on this node as well. */
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_gradient);
+ scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_transport);
+#endif
}
/* Is the foreign cell active and will need stuff from us? */
@@ -1430,6 +1487,16 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
scheduler_activate_send(s, ci->mpi.send,
task_subtype_rt_transport, cj_nodeID);
}
+ } else if (ci_active_rt) {
+#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
+ /* If the foreign cell is inactive, but the local cell is active,
+ * we still need to send stuff to be able to do the force
+ * interaction on both nodes */
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_rt_gradient,
+ cj_nodeID);
+ scheduler_activate_send(s, ci->mpi.send, task_subtype_rt_transport,
+ cj_nodeID);
+#endif
}
}
#endif
diff --git a/src/hydro/Gizmo/MFM/hydro_part.h b/src/hydro/Gizmo/MFM/hydro_part.h
index 509dd79ad661317714425bfc93939264f68c5b33..f476b860eaaca8d54386fe3f2502db3b58c67df5 100644
--- a/src/hydro/Gizmo/MFM/hydro_part.h
+++ b/src/hydro/Gizmo/MFM/hydro_part.h
@@ -84,19 +84,6 @@ struct part {
} limiter;
struct {
- /* Variables used for timestep calculation. */
- struct {
-
- /* Maximum signal velocity among all the neighbours of the particle. The
- * signal velocity encodes information about the relative fluid
- * velocities
- * AND particle velocities of the neighbour and this particle, as well
- * as
- * the sound speed of both particles. */
- float vmax;
-
- } timestepvars;
-
/* Quantities used during the force loop. */
struct {
@@ -165,6 +152,19 @@ struct part {
} geometry;
+ /* Variables used for timestep calculation. */
+ struct {
+
+ /* Maximum signal velocity among all the neighbours of the particle. The
+ * signal velocity encodes information about the relative fluid
+ * velocities
+ * AND particle velocities of the neighbour and this particle, as well
+ * as
+ * the sound speed of both particles. */
+ float vmax;
+
+ } timestepvars;
+
/*! Chemistry information */
struct chemistry_part_data chemistry_data;
diff --git a/src/hydro/Gizmo/hydro_getters.h b/src/hydro/Gizmo/hydro_getters.h
index c97d20d1b3855fe8f98565822df3fb945ae261b9..7a83e093680ed63466eedd8c3885cf19a9279b5c 100644
--- a/src/hydro/Gizmo/hydro_getters.h
+++ b/src/hydro/Gizmo/hydro_getters.h
@@ -113,7 +113,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_get_slope_limiter(
* @param p The particle of interest.
*/
__attribute__((always_inline)) INLINE static float
-hydro_get_comoving_internal_energy(const struct part* restrict p) {
+hydro_get_comoving_internal_energy(const struct part* restrict p,
+ const struct xpart* restrict xp) {
if (p->rho > 0.0f)
return gas_internal_energy_from_pressure(p->rho, p->P);
@@ -134,7 +135,7 @@ hydro_get_physical_internal_energy(const struct part* restrict p,
const struct cosmology* cosmo) {
return cosmo->a_factor_internal_energy *
- hydro_get_comoving_internal_energy(p);
+ hydro_get_comoving_internal_energy(p, xp);
}
/**
@@ -158,7 +159,7 @@ hydro_get_drifted_physical_internal_energy(const struct part* restrict p,
__attribute__((always_inline)) INLINE static float
hydro_get_drifted_comoving_internal_energy(const struct part* restrict p) {
- return hydro_get_comoving_internal_energy(p);
+ return hydro_get_comoving_internal_energy(p, NULL);
}
/**
@@ -322,8 +323,30 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_internal_energy_dt(const struct part* restrict p) {
- error("Needs implementing");
- return 0.0f;
+ float W[5];
+ hydro_part_get_primitive_variables(p, W);
+
+ if (W[0] <= 0.0f) {
+ return 0.0f;
+ }
+
+ const float rho_inv = 1.f / W[0];
+
+ float gradrho[3], gradvx[3], gradvy[3], gradvz[3], gradP[3];
+ hydro_part_get_gradients(p, gradrho, gradvx, gradvy, gradvz, gradP);
+
+ const float divv = gradvx[0] + gradvy[1] + gradvz[2];
+
+ float gradu[3] = {0.f, 0.f, 0.f};
+ for (int i = 0; i < 3; i++) {
+ gradu[i] = hydro_one_over_gamma_minus_one * rho_inv *
+ (gradP[i] - rho_inv * W[4] * gradrho[i]);
+ }
+
+ const float du_dt = -(W[1] * gradu[0] + W[2] * gradu[1] + W[3] * gradu[2]) -
+ rho_inv * W[4] * divv;
+
+ return du_dt;
}
/**
@@ -337,8 +360,9 @@ hydro_get_comoving_internal_energy_dt(const struct part* restrict p) {
__attribute__((always_inline)) INLINE static float
hydro_get_physical_internal_energy_dt(const struct part* restrict p,
const struct cosmology* cosmo) {
- error("Needs implementing");
- return 0.0f;
+
+ return hydro_get_comoving_internal_energy_dt(p) *
+ cosmo->a_factor_internal_energy;
}
/**
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index c56db563e02d31cfdf8b99b99a9ade21d65c1796..171484205512e86718f4c5de60d826b0c77c7985 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -76,7 +76,7 @@ INLINE static void hydro_read_particles(struct part* parts,
INLINE static void convert_u(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
- ret[0] = hydro_get_comoving_internal_energy(p);
+ ret[0] = hydro_get_comoving_internal_energy(p, xp);
}
/**
diff --git a/src/hydro/Gizmo/hydro_setters.h b/src/hydro/Gizmo/hydro_setters.h
index fd7e6f8ca35cd0879bf2b0e2f40eba41fbb257d3..640912ccbb226a74f629afb0ac630079f6332854 100644
--- a/src/hydro/Gizmo/hydro_setters.h
+++ b/src/hydro/Gizmo/hydro_setters.h
@@ -195,7 +195,10 @@ __attribute__((always_inline)) INLINE static void hydro_set_mass(
__attribute__((always_inline)) INLINE static void
hydro_set_comoving_internal_energy_dt(struct part* restrict p,
const float du_dt) {
- error("Needs implementing");
+
+ const float old_du_dt = hydro_get_comoving_internal_energy_dt(p);
+
+ p->flux.energy += p->conserved.mass * (du_dt - old_du_dt) * p->flux.dt;
}
/**
@@ -211,8 +214,43 @@ __attribute__((always_inline)) INLINE static void
hydro_set_physical_internal_energy_dt(struct part* restrict p,
const struct cosmology* restrict cosmo,
const float du_dt) {
- error("Needs implementing");
+
+ hydro_set_comoving_internal_energy_dt(
+ p, du_dt / cosmo->a_factor_internal_energy);
}
+
+/**
+ * @brief Sets the comoving internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param u The comoving internal energy
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_comoving_internal_energy(struct part* p, const float u) {
+
+ const float mass = p->conserved.mass;
+ if (mass <= 0.0f) {
+ return;
+ }
+
+ const float Etherm = mass * u;
+
+#ifdef GIZMO_TOTAL_ENERGY
+ const float Ekin = 0.5f *
+ (p->conserved.momentum[0] * p->conserved.momentum[0] +
+ p->conserved.momentum[1] * p->conserved.momentum[1] +
+ p->conserved.momentum[2] * p->conserved.momentum[2]) /
+ mass;
+
+ const float Etot = Ekin + Etherm;
+ p->conserved.energy = Etot;
+#else
+ p->conserved.energy = Etherm;
+#endif
+
+ p->P = gas_pressure_from_internal_energy(p->rho, u);
+}
+
/**
* @brief Sets the physical entropy of a particle
*
@@ -225,7 +263,8 @@ __attribute__((always_inline)) INLINE static void hydro_set_physical_entropy(
struct part* p, struct xpart* xp, const struct cosmology* cosmo,
const float entropy) {
- error("Needs implementing");
+ const float u = gas_internal_energy_from_entropy(p->rho, entropy);
+ hydro_set_comoving_internal_energy(p, u);
}
/**
@@ -240,7 +279,8 @@ __attribute__((always_inline)) INLINE static void
hydro_set_physical_internal_energy(struct part* p, struct xpart* xp,
const struct cosmology* cosmo,
const float u) {
- error("Need implementing");
+
+ hydro_set_comoving_internal_energy(p, u / cosmo->a_factor_internal_energy);
}
/**
@@ -254,7 +294,8 @@ __attribute__((always_inline)) INLINE static void
hydro_set_drifted_physical_internal_energy(struct part* p,
const struct cosmology* cosmo,
const float u) {
- error("Need implementing");
+
+ hydro_set_comoving_internal_energy(p, u / cosmo->a_factor_internal_energy);
}
/**
diff --git a/src/multipole.h b/src/multipole.h
index 897623f9515add2e40f0cc544b1e9702ddddf96f..26c163a7c952361c63f923a60ca040647203a57e 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1031,7 +1031,7 @@ __attribute__((nonnull)) INLINE static void gravity_P2M(
float min_delta_vel[3] = {0., 0., 0.};
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
- double M_100 = 0., M_010 = 0., M_001 = 0.;
+ /* double M_100 = 0., M_010 = 0., M_001 = 0.; */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
double M_200 = 0., M_020 = 0., M_002 = 0.;
@@ -1086,9 +1086,9 @@ __attribute__((nonnull)) INLINE static void gravity_P2M(
const double m = gparts[k].mass;
/* 1st order terms */
- M_100 += -m * X_100(dx);
- M_010 += -m * X_010(dx);
- M_001 += -m * X_001(dx);
+ /* M_100 += -m * X_100(dx); */
+ /* M_010 += -m * X_010(dx); */
+ /* M_001 += -m * X_001(dx); */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
diff --git a/src/power_spectrum.c b/src/power_spectrum.c
index 0b947fee9e65ecb6f55ac044b526bd1bde80f2e7..324d48c231b6a6189973a40476640e4ec62461a2 100644
--- a/src/power_spectrum.c
+++ b/src/power_spectrum.c
@@ -51,6 +51,8 @@
#define power_data_default_fold_factor 4
#define power_data_default_window_order 3
+#ifdef HAVE_FFTW
+
/**
* @brief Return the #power_type corresponding to a given string.
*/
@@ -104,137 +106,6 @@ INLINE static const char* get_powtype_filename(const enum power_type type) {
return powtype_filenames[type];
}
-/**
- * @brief Initialize power spectra calculation.
- *
- * Reads in the power spectrum parameters, sets up FFTW
- * (likely already done for PM), then sets up an FFT plan
- * that can be reused every time.
- *
- * @param nr_threads The number of threads used.
- */
-void power_init(struct power_spectrum_data* p, struct swift_params* params,
- int nr_threads) {
-
-#ifdef HAVE_FFTW
-
- /* Get power spectrum parameters */
- p->Ngrid = parser_get_opt_param_int(params, "PowerSpectrum:grid_side_length",
- power_data_default_grid_side_length);
- p->Nfold = parser_get_param_int(params, "PowerSpectrum:num_folds");
- p->foldfac = parser_get_opt_param_int(params, "PowerSpectrum:fold_factor",
- power_data_default_fold_factor);
- p->windoworder = parser_get_opt_param_int(
- params, "PowerSpectrum:window_order", power_data_default_window_order);
-
- if (p->windoworder > 3 || p->windoworder < 1)
- error("Power spectrum calculation is not implemented for %dth order!",
- p->windoworder);
- if (p->windoworder == 1)
- message("WARNING: window order is recommended to be at least 2 (CIC).");
- if (p->windoworder <= 2 && p->foldfac > 4)
- message(
- "WARNING: fold factor is recommended not to exceed 4 for a "
- "mass assignment order of 2 (CIC) or below.");
- if (p->windoworder == 3 && p->foldfac > 6)
- message(
- "WARNING: fold factor is recommended not to exceed 6 for a "
- "mass assignment order of 3 (TSC) or below.");
-
- /* Make sensible choices for the k-cuts */
- const int kcutn = (p->windoworder >= 3) ? 90 : 70;
- const int kcutleft = (int)(p->Ngrid / 256.0 * kcutn);
- const int kcutright = (int)(p->Ngrid / 256.0 * (double)kcutn / p->foldfac);
- if (kcutright < 10 || (kcutleft - kcutright) < 30)
- error(
- "Combination of power grid size and fold factor do not allow for "
- "enough overlap between foldings!");
-
- p->nr_threads = nr_threads;
-
-#ifdef HAVE_THREADED_FFTW
- /* Initialise the thread-parallel FFTW version
- (probably already done for the PM, but does not matter) */
- if (p->Ngrid >= 64) {
- fftw_init_threads();
- fftw_plan_with_nthreads(nr_threads);
- }
-#else
- message("Note that FFTW is not threaded!");
-#endif
-
- char** requested_spectra;
- parser_get_param_string_array(params, "PowerSpectrum:requested_spectra",
- &p->spectrumcount, &requested_spectra);
-
- p->types1 =
- (enum power_type*)malloc(p->spectrumcount * sizeof(enum power_type));
- p->types2 =
- (enum power_type*)malloc(p->spectrumcount * sizeof(enum power_type));
-
- /* Parse which spectra are being requested */
- for (int i = 0; i < p->spectrumcount; ++i) {
-
- char* pstr = strtok(requested_spectra[i], "-");
- if (pstr == NULL)
- error("Requested power spectra are not in the format type1-type2!");
- char type1[32];
- strcpy(type1, pstr);
-
- pstr = strtok(NULL, "-");
- if (pstr == NULL)
- error("Requested power spectra are not in the format type1-type2!");
- char type2[32];
- strcpy(type2, pstr);
-
- p->types1[i] = power_spectrum_get_type(type1);
- p->types2[i] = power_spectrum_get_type(type2);
- }
-
- /* Initialize the plan only once -- much faster for FFTs run often!
- * Does require us to allocate the grids, but we delete them right away.
- * Plan can only be used for the same FFTW call */
- const int Ngrid = p->Ngrid;
-
- /* Grid is padded to allow for in-place FFT */
- p->powgrid = fftw_alloc_real(Ngrid * Ngrid * (Ngrid + 2));
- /* Pointer to grid to interpret it as complex data */
- p->powgridft = (fftw_complex*)p->powgrid;
-
- p->fftplanpow = fftw_plan_dft_r2c_3d(Ngrid, Ngrid, Ngrid, p->powgrid,
- p->powgridft, FFTW_MEASURE);
-
- fftw_free(p->powgrid);
- p->powgrid = NULL;
- p->powgridft = NULL;
-
- /* Do the same for a second grid/plan to allow for cross power */
-
- /* Grid is padded to allow for in-place FFT */
- p->powgrid2 = fftw_alloc_real(Ngrid * Ngrid * (Ngrid + 2));
- /* Pointer to grid to interpret it as complex data */
- p->powgridft2 = (fftw_complex*)p->powgrid2;
-
- p->fftplanpow2 = fftw_plan_dft_r2c_3d(Ngrid, Ngrid, Ngrid, p->powgrid2,
- p->powgridft2, FFTW_MEASURE);
-
- fftw_free(p->powgrid2);
- p->powgrid2 = NULL;
- p->powgridft2 = NULL;
-
- /* Create directories for power spectra and foldings */
- if (engine_rank == 0) {
- safe_checkdir("power_spectra", /*create=*/1);
- safe_checkdir("power_spectra/foldings", /*create=*/1);
- }
-
-#else
- error("Trying to initialize the PS code without FFTW present!");
-#endif
-}
-
-#ifdef HAVE_FFTW
-
/**
* @brief Shared information for shot noise to be used by all the threads in the
* pool.
@@ -1365,6 +1236,135 @@ void power_spectrum(const enum power_type type1, const enum power_type type2,
#endif /* HAVE_FFTW */
+/**
+ * @brief Initialize power spectra calculation.
+ *
+ * Reads in the power spectrum parameters, sets up FFTW
+ * (likely already done for PM), then sets up an FFT plan
+ * that can be reused every time.
+ *
+ * @param nr_threads The number of threads used.
+ */
+void power_init(struct power_spectrum_data* p, struct swift_params* params,
+ int nr_threads) {
+
+#ifdef HAVE_FFTW
+
+ /* Get power spectrum parameters */
+ p->Ngrid = parser_get_opt_param_int(params, "PowerSpectrum:grid_side_length",
+ power_data_default_grid_side_length);
+ p->Nfold = parser_get_param_int(params, "PowerSpectrum:num_folds");
+ p->foldfac = parser_get_opt_param_int(params, "PowerSpectrum:fold_factor",
+ power_data_default_fold_factor);
+ p->windoworder = parser_get_opt_param_int(
+ params, "PowerSpectrum:window_order", power_data_default_window_order);
+
+ if (p->windoworder > 3 || p->windoworder < 1)
+ error("Power spectrum calculation is not implemented for %dth order!",
+ p->windoworder);
+ if (p->windoworder == 1)
+ message("WARNING: window order is recommended to be at least 2 (CIC).");
+ if (p->windoworder <= 2 && p->foldfac > 4)
+ message(
+ "WARNING: fold factor is recommended not to exceed 4 for a "
+ "mass assignment order of 2 (CIC) or below.");
+ if (p->windoworder == 3 && p->foldfac > 6)
+ message(
+ "WARNING: fold factor is recommended not to exceed 6 for a "
+ "mass assignment order of 3 (TSC) or below.");
+
+ /* Make sensible choices for the k-cuts */
+ const int kcutn = (p->windoworder >= 3) ? 90 : 70;
+ const int kcutleft = (int)(p->Ngrid / 256.0 * kcutn);
+ const int kcutright = (int)(p->Ngrid / 256.0 * (double)kcutn / p->foldfac);
+ if (kcutright < 10 || (kcutleft - kcutright) < 30)
+ error(
+ "Combination of power grid size and fold factor do not allow for "
+ "enough overlap between foldings!");
+
+ p->nr_threads = nr_threads;
+
+#ifdef HAVE_THREADED_FFTW
+ /* Initialise the thread-parallel FFTW version
+ (probably already done for the PM, but does not matter) */
+ if (p->Ngrid >= 64) {
+ fftw_init_threads();
+ fftw_plan_with_nthreads(nr_threads);
+ }
+#else
+ message("Note that FFTW is not threaded!");
+#endif
+
+ char** requested_spectra;
+ parser_get_param_string_array(params, "PowerSpectrum:requested_spectra",
+ &p->spectrumcount, &requested_spectra);
+
+ p->types1 =
+ (enum power_type*)malloc(p->spectrumcount * sizeof(enum power_type));
+ p->types2 =
+ (enum power_type*)malloc(p->spectrumcount * sizeof(enum power_type));
+
+ /* Parse which spectra are being requested */
+ for (int i = 0; i < p->spectrumcount; ++i) {
+
+ char* pstr = strtok(requested_spectra[i], "-");
+ if (pstr == NULL)
+ error("Requested power spectra are not in the format type1-type2!");
+ char type1[32];
+ strcpy(type1, pstr);
+
+ pstr = strtok(NULL, "-");
+ if (pstr == NULL)
+ error("Requested power spectra are not in the format type1-type2!");
+ char type2[32];
+ strcpy(type2, pstr);
+
+ p->types1[i] = power_spectrum_get_type(type1);
+ p->types2[i] = power_spectrum_get_type(type2);
+ }
+
+ /* Initialize the plan only once -- much faster for FFTs run often!
+ * Does require us to allocate the grids, but we delete them right away.
+ * Plan can only be used for the same FFTW call */
+ const int Ngrid = p->Ngrid;
+
+ /* Grid is padded to allow for in-place FFT */
+ p->powgrid = fftw_alloc_real(Ngrid * Ngrid * (Ngrid + 2));
+ /* Pointer to grid to interpret it as complex data */
+ p->powgridft = (fftw_complex*)p->powgrid;
+
+ p->fftplanpow = fftw_plan_dft_r2c_3d(Ngrid, Ngrid, Ngrid, p->powgrid,
+ p->powgridft, FFTW_MEASURE);
+
+ fftw_free(p->powgrid);
+ p->powgrid = NULL;
+ p->powgridft = NULL;
+
+ /* Do the same for a second grid/plan to allow for cross power */
+
+ /* Grid is padded to allow for in-place FFT */
+ p->powgrid2 = fftw_alloc_real(Ngrid * Ngrid * (Ngrid + 2));
+ /* Pointer to grid to interpret it as complex data */
+ p->powgridft2 = (fftw_complex*)p->powgrid2;
+
+ p->fftplanpow2 = fftw_plan_dft_r2c_3d(Ngrid, Ngrid, Ngrid, p->powgrid2,
+ p->powgridft2, FFTW_MEASURE);
+
+ fftw_free(p->powgrid2);
+ p->powgrid2 = NULL;
+ p->powgridft2 = NULL;
+
+ /* Create directories for power spectra and foldings */
+ if (engine_rank == 0) {
+ safe_checkdir("power_spectra", /*create=*/1);
+ safe_checkdir("power_spectra/foldings", /*create=*/1);
+ }
+
+#else
+ error("Trying to initialize the PS code without FFTW present!");
+#endif
+}
+
void calc_all_power_spectra(struct power_spectrum_data* pow_data,
const struct space* s, struct threadpool* tp,
const int verbose) {
diff --git a/src/rt/GEAR/rt.h b/src/rt/GEAR/rt.h
index 92d8c3d11e954033963d805c573de9b64c939b20..58e2e8f0ec6800167b6a74413a5e60966cb73d77 100644
--- a/src/rt/GEAR/rt.h
+++ b/src/rt/GEAR/rt.h
@@ -132,9 +132,11 @@ __attribute__((always_inline)) INLINE static void rt_reset_part(
* @brief Reset RT particle data which needs to be reset each sub-cycle.
*
* @param p the particle to work on
+ * @param cosmo Cosmology.
+ * @param dt the current particle RT time step
*/
__attribute__((always_inline)) INLINE static void rt_reset_part_each_subcycle(
- struct part* restrict p) {
+ struct part* restrict p, const struct cosmology* cosmo, double dt) {
#ifdef SWIFT_RT_DEBUG_CHECKS
rt_debugging_reset_each_subcycle(p);
@@ -144,7 +146,8 @@ __attribute__((always_inline)) INLINE static void rt_reset_part_each_subcycle(
/* the Gizmo-style slope limiting doesn't help for RT as is,
* so we're skipping it for now. */
/* rt_slope_limit_cell_init(p); */
- rt_part_reset_fluxes(p);
+
+ p->rt_data.flux_dt = dt;
}
/**
@@ -162,7 +165,8 @@ __attribute__((always_inline)) INLINE static void rt_first_init_part(
rt_init_part(p);
rt_reset_part(p, cosmo);
rt_part_reset_mass_fluxes(p);
- rt_reset_part_each_subcycle(p);
+ rt_reset_part_each_subcycle(p, cosmo, 0.);
+ rt_part_reset_fluxes(p);
#ifdef SWIFT_RT_DEBUG_CHECKS
p->rt_data.debug_radiation_absorbed_tot = 0ULL;
@@ -477,15 +481,20 @@ __attribute__((always_inline)) INLINE static void rt_finalise_transport(
for (int g = 0; g < RT_NGROUPS; g++) {
const float e_old = rtd->radiation[g].energy_density;
- /* Note: in this scheme, we're updating d/dt (U * V) + sum F * A = 0.
+ /* Note: in this scheme, we're updating d/dt (U * V) + sum F * A * dt = 0.
* So we'll need the division by the volume here. */
- rtd->radiation[g].energy_density += rtd->flux[g].energy * dt * Vinv;
- rtd->radiation[g].flux[0] += rtd->flux[g].flux[0] * dt * Vinv;
- rtd->radiation[g].flux[1] += rtd->flux[g].flux[1] * dt * Vinv;
- rtd->radiation[g].flux[2] += rtd->flux[g].flux[2] * dt * Vinv;
+ rtd->radiation[g].energy_density += rtd->flux[g].energy * Vinv;
+ rtd->radiation[g].flux[0] += rtd->flux[g].flux[0] * Vinv;
+ rtd->radiation[g].flux[1] += rtd->flux[g].flux[1] * Vinv;
+ rtd->radiation[g].flux[2] += rtd->flux[g].flux[2] * Vinv;
rt_check_unphysical_state(&rtd->radiation[g].energy_density,
rtd->radiation[g].flux, e_old, /*callloc=*/4);
}
+
+ /* Reset the fluxes now that they have been applied. */
+ rt_part_reset_fluxes(p);
+ /* Mark the particle as inactive for now. */
+ rtd->flux_dt = -1.f;
}
/**
@@ -549,54 +558,60 @@ __attribute__((always_inline)) INLINE static void rt_kick_extra(
}
#endif
- /* Update the mass fraction changes due to interparticle fluxes */
- const float current_mass = p->conserved.mass;
-
- const float current_mass_HI =
- current_mass * p->rt_data.tchem.mass_fraction_HI;
- const float current_mass_HII =
- current_mass * p->rt_data.tchem.mass_fraction_HII;
- const float current_mass_HeI =
- current_mass * p->rt_data.tchem.mass_fraction_HeI;
- const float current_mass_HeII =
- current_mass * p->rt_data.tchem.mass_fraction_HeII;
- const float current_mass_HeIII =
- current_mass * p->rt_data.tchem.mass_fraction_HeIII;
-
- const float new_mass_HI =
- current_mass_HI + p->rt_data.mass_flux.HI * dt_therm;
- const float new_mass_HII =
- current_mass_HII + p->rt_data.mass_flux.HII * dt_therm;
- const float new_mass_HeI =
- current_mass_HeI + p->rt_data.mass_flux.HeI * dt_therm;
- const float new_mass_HeII =
- current_mass_HeII + p->rt_data.mass_flux.HeII * dt_therm;
- const float new_mass_HeIII =
- current_mass_HeIII + p->rt_data.mass_flux.HeIII * dt_therm;
-
- const float new_mass_tot = new_mass_HI + new_mass_HII + new_mass_HeI +
- new_mass_HeII + new_mass_HeIII;
-
- /* During the initial fake time step, if the mass fractions haven't been set
- * up yet, we could encounter divisions by zero here, so skip that. If it
- * isn't the initial time step, the error will be caught down the line by
- * another call
- * to rt_check_unphysical_mass_fractions() (not the one 10 lines below this)
- */
- if (new_mass_tot == 0.f) return;
-
- const float new_mass_tot_inv = 1.f / new_mass_tot;
-
- p->rt_data.tchem.mass_fraction_HI = new_mass_HI * new_mass_tot_inv;
- p->rt_data.tchem.mass_fraction_HII = new_mass_HII * new_mass_tot_inv;
- p->rt_data.tchem.mass_fraction_HeI = new_mass_HeI * new_mass_tot_inv;
- p->rt_data.tchem.mass_fraction_HeII = new_mass_HeII * new_mass_tot_inv;
- p->rt_data.tchem.mass_fraction_HeIII = new_mass_HeIII * new_mass_tot_inv;
-
- rt_check_unphysical_mass_fractions(p);
-
- /* Don't update actual particle mass, that'll be done in the
- * hydro_kick_extra calls */
+ /* Note: We need to mimick here what Gizmo does for the mass fluxes.
+ * The relevant time scale is the hydro time step for the mass fluxes,
+ * not the RT times. We also need to prevent the kick to apply the mass
+ * fluxes twice, so exit if the particle time step < 0 */
+
+ if (p->flux.dt > 0.0f) {
+ /* Update the mass fraction changes due to interparticle fluxes */
+ const float current_mass = p->conserved.mass;
+
+ const float current_mass_HI =
+ current_mass * p->rt_data.tchem.mass_fraction_HI;
+ const float current_mass_HII =
+ current_mass * p->rt_data.tchem.mass_fraction_HII;
+ const float current_mass_HeI =
+ current_mass * p->rt_data.tchem.mass_fraction_HeI;
+ const float current_mass_HeII =
+ current_mass * p->rt_data.tchem.mass_fraction_HeII;
+ const float current_mass_HeIII =
+ current_mass * p->rt_data.tchem.mass_fraction_HeIII;
+
+ const float new_mass_HI = current_mass_HI + p->rt_data.mass_flux.HI;
+ const float new_mass_HII = current_mass_HII + p->rt_data.mass_flux.HII;
+ const float new_mass_HeI = current_mass_HeI + p->rt_data.mass_flux.HeI;
+ const float new_mass_HeII = current_mass_HeII + p->rt_data.mass_flux.HeII;
+ const float new_mass_HeIII =
+ current_mass_HeIII + p->rt_data.mass_flux.HeIII;
+
+ const float new_mass_tot = new_mass_HI + new_mass_HII + new_mass_HeI +
+ new_mass_HeII + new_mass_HeIII;
+
+ /* During the initial fake time step, if the mass fractions haven't been set
+ * up yet, we could encounter divisions by zero here, so skip that. If it
+ * isn't the initial time step, the error will be caught down the line by
+ * another call to rt_check_unphysical_mass_fractions() (not the one 10
+ * lines below this) */
+ if (new_mass_tot == 0.f) return;
+
+ const float new_mass_tot_inv = 1.f / new_mass_tot;
+
+ p->rt_data.tchem.mass_fraction_HI = new_mass_HI * new_mass_tot_inv;
+ p->rt_data.tchem.mass_fraction_HII = new_mass_HII * new_mass_tot_inv;
+ p->rt_data.tchem.mass_fraction_HeI = new_mass_HeI * new_mass_tot_inv;
+ p->rt_data.tchem.mass_fraction_HeII = new_mass_HeII * new_mass_tot_inv;
+ p->rt_data.tchem.mass_fraction_HeIII = new_mass_HeIII * new_mass_tot_inv;
+
+ rt_check_unphysical_mass_fractions(p);
+
+ /* Reset fluxes after they have been applied, so they can be collected
+ * again even when particle is inactive. */
+ rt_part_reset_mass_fluxes(p);
+
+ /* Don't update actual particle mass, that'll be done in the
+ * hydro_kick_extra calls */
+ }
}
/**
@@ -609,10 +624,7 @@ __attribute__((always_inline)) INLINE static void rt_kick_extra(
* @param p particle to work on
**/
__attribute__((always_inline)) INLINE static void rt_prepare_force(
- struct part* p) {
-
- rt_part_reset_mass_fluxes(p);
-}
+ struct part* p) {}
/**
* @brief Extra operations to be done during the drift
diff --git a/src/rt/GEAR/rt_additions.h b/src/rt/GEAR/rt_additions.h
index 34baa6ac94797afa1099e72d291af986aed85b1d..6e2e269c3cca08b23d547c992006d9b3407c8b51 100644
--- a/src/rt/GEAR/rt_additions.h
+++ b/src/rt/GEAR/rt_additions.h
@@ -43,58 +43,49 @@ __attribute__((always_inline)) INLINE static void rt_part_update_mass_fluxes(
* its own mass fractions. If it's gaining mass, it's gaining mass
* according to the interacting particle's mass fractions. */
+ /* get the time step for the flux exchange. This is always the smallest time
+ step among the two particles */
+ const float mindt =
+ (pj->flux.dt > 0.0f) ? fminf(pi->flux.dt, pj->flux.dt) : pi->flux.dt;
+
+ const float mdt = mass_flux * mindt;
+
/* Convention: negative flux for "left" particle pi */
if (mass_flux < 0.f) {
/* "left" particle is gaining mass */
- pi->rt_data.mass_flux.HI -= pj->rt_data.tchem.mass_fraction_HI * mass_flux;
- pi->rt_data.mass_flux.HII -=
- pj->rt_data.tchem.mass_fraction_HII * mass_flux;
- pi->rt_data.mass_flux.HeI -=
- pj->rt_data.tchem.mass_fraction_HeI * mass_flux;
- pi->rt_data.mass_flux.HeII -=
- pj->rt_data.tchem.mass_fraction_HeII * mass_flux;
- pi->rt_data.mass_flux.HeIII -=
- pj->rt_data.tchem.mass_fraction_HeIII * mass_flux;
+ pi->rt_data.mass_flux.HI -= pj->rt_data.tchem.mass_fraction_HI * mdt;
+ pi->rt_data.mass_flux.HII -= pj->rt_data.tchem.mass_fraction_HII * mdt;
+ pi->rt_data.mass_flux.HeI -= pj->rt_data.tchem.mass_fraction_HeI * mdt;
+ pi->rt_data.mass_flux.HeII -= pj->rt_data.tchem.mass_fraction_HeII * mdt;
+ pi->rt_data.mass_flux.HeIII -= pj->rt_data.tchem.mass_fraction_HeIII * mdt;
} else {
/* "left" particle is losing mass */
- pi->rt_data.mass_flux.HI -= pi->rt_data.tchem.mass_fraction_HI * mass_flux;
- pi->rt_data.mass_flux.HII -=
- pi->rt_data.tchem.mass_fraction_HII * mass_flux;
- pi->rt_data.mass_flux.HeI -=
- pi->rt_data.tchem.mass_fraction_HeI * mass_flux;
- pi->rt_data.mass_flux.HeII -=
- pi->rt_data.tchem.mass_fraction_HeII * mass_flux;
- pi->rt_data.mass_flux.HeIII -=
- pi->rt_data.tchem.mass_fraction_HeIII * mass_flux;
+ pi->rt_data.mass_flux.HI -= pi->rt_data.tchem.mass_fraction_HI * mdt;
+ pi->rt_data.mass_flux.HII -= pi->rt_data.tchem.mass_fraction_HII * mdt;
+ pi->rt_data.mass_flux.HeI -= pi->rt_data.tchem.mass_fraction_HeI * mdt;
+ pi->rt_data.mass_flux.HeII -= pi->rt_data.tchem.mass_fraction_HeII * mdt;
+ pi->rt_data.mass_flux.HeIII -= pi->rt_data.tchem.mass_fraction_HeIII * mdt;
}
- if (mode == 1) {
- /* Update right particle as well */
+ if (mode == 1 || (pj->flux.dt < 0.f)) {
+ /* Update right particle as well, even if it's inactive */
if (mass_flux > 0.f) {
/* "right" particle is gaining mass */
- pj->rt_data.mass_flux.HI +=
- pi->rt_data.tchem.mass_fraction_HI * mass_flux;
- pj->rt_data.mass_flux.HII +=
- pi->rt_data.tchem.mass_fraction_HII * mass_flux;
- pj->rt_data.mass_flux.HeI +=
- pi->rt_data.tchem.mass_fraction_HeI * mass_flux;
- pj->rt_data.mass_flux.HeII +=
- pi->rt_data.tchem.mass_fraction_HeII * mass_flux;
+ pj->rt_data.mass_flux.HI += pi->rt_data.tchem.mass_fraction_HI * mdt;
+ pj->rt_data.mass_flux.HII += pi->rt_data.tchem.mass_fraction_HII * mdt;
+ pj->rt_data.mass_flux.HeI += pi->rt_data.tchem.mass_fraction_HeI * mdt;
+ pj->rt_data.mass_flux.HeII += pi->rt_data.tchem.mass_fraction_HeII * mdt;
pj->rt_data.mass_flux.HeIII +=
- pi->rt_data.tchem.mass_fraction_HeIII * mass_flux;
+ pi->rt_data.tchem.mass_fraction_HeIII * mdt;
} else {
/* "right" particle is losing mass */
- pj->rt_data.mass_flux.HI +=
- pj->rt_data.tchem.mass_fraction_HI * mass_flux;
- pj->rt_data.mass_flux.HII +=
- pj->rt_data.tchem.mass_fraction_HII * mass_flux;
- pj->rt_data.mass_flux.HeI +=
- pj->rt_data.tchem.mass_fraction_HeI * mass_flux;
- pj->rt_data.mass_flux.HeII +=
- pj->rt_data.tchem.mass_fraction_HeII * mass_flux;
+ pj->rt_data.mass_flux.HI += pj->rt_data.tchem.mass_fraction_HI * mdt;
+ pj->rt_data.mass_flux.HII += pj->rt_data.tchem.mass_fraction_HII * mdt;
+ pj->rt_data.mass_flux.HeI += pj->rt_data.tchem.mass_fraction_HeI * mdt;
+ pj->rt_data.mass_flux.HeII += pj->rt_data.tchem.mass_fraction_HeII * mdt;
pj->rt_data.mass_flux.HeIII +=
- pj->rt_data.tchem.mass_fraction_HeIII * mass_flux;
+ pj->rt_data.tchem.mass_fraction_HeIII * mdt;
}
}
}
diff --git a/src/rt/GEAR/rt_iact.h b/src/rt/GEAR/rt_iact.h
index 927d5881f9e7fb0a2d23fbe16e4532afe9af3c67..3ad71d3b66616369ee45851cca73b4cbb9fae0ed 100644
--- a/src/rt/GEAR/rt_iact.h
+++ b/src/rt/GEAR/rt_iact.h
@@ -326,6 +326,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_rt_flux_common(
struct rt_part_data *restrict rti = &pi->rt_data;
struct rt_part_data *restrict rtj = &pj->rt_data;
+ /* Get the time step for the flux exchange. This is always the smallest time
+ * step among the two particles. */
+ const float mindt =
+ (rtj->flux_dt > 0.f) ? fminf(rti->flux_dt, rtj->flux_dt) : rti->flux_dt;
for (int g = 0; g < RT_NGROUPS; g++) {
@@ -350,15 +354,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_rt_flux_common(
* flux is subtracted from the left state, and added to the right
* state, based on how we chose the unit vector. By this convention,
* the time integration results in conserved quantity += flux * dt */
- rti->flux[g].energy -= totflux[0];
- rti->flux[g].flux[0] -= totflux[1];
- rti->flux[g].flux[1] -= totflux[2];
- rti->flux[g].flux[2] -= totflux[3];
- if (mode == 1) {
- rtj->flux[g].energy += totflux[0];
- rtj->flux[g].flux[0] += totflux[1];
- rtj->flux[g].flux[1] += totflux[2];
- rtj->flux[g].flux[2] += totflux[3];
+ /* Unlike in SPH schemes, we do need to update inactive neighbours, so that
+ * the fluxes are always exchanged symmetrically. Thanks to our sneaky use
+ * of flux_dt, we can detect inactive neighbours through their negative time
+ * step. */
+ rti->flux[g].energy -= totflux[0] * mindt;
+ rti->flux[g].flux[0] -= totflux[1] * mindt;
+ rti->flux[g].flux[1] -= totflux[2] * mindt;
+ rti->flux[g].flux[2] -= totflux[3] * mindt;
+ if (mode == 1 || (rtj->flux_dt < 0.f)) {
+ rtj->flux[g].energy += totflux[0] * mindt;
+ rtj->flux[g].flux[0] += totflux[1] * mindt;
+ rtj->flux[g].flux[1] += totflux[2] * mindt;
+ rtj->flux[g].flux[2] += totflux[3] * mindt;
}
}
}
diff --git a/src/rt/GEAR/rt_struct.h b/src/rt/GEAR/rt_struct.h
index 632bc15b759588007682c6c71b67bdbe18221f65..63cfebbb6a2e3337e317edea206d7dd183f7ac6e 100644
--- a/src/rt/GEAR/rt_struct.h
+++ b/src/rt/GEAR/rt_struct.h
@@ -39,6 +39,9 @@ struct rt_part_data {
float flux[3];
} flux[RT_NGROUPS];
+ /* Particle RT time step. */
+ float flux_dt;
+
/* gradients of the radiation state. */
/* for the flux[3][3] quantity:
* first index: x, y, z coordinate of the flux.
diff --git a/src/rt/SPHM1RT/rt.h b/src/rt/SPHM1RT/rt.h
index 0ae7c5805284fb83e684dbe75d5a3c373a7883f5..88437dfc1fe657f5e51037935fbf527484ffe596 100644
--- a/src/rt/SPHM1RT/rt.h
+++ b/src/rt/SPHM1RT/rt.h
@@ -55,7 +55,17 @@ __attribute__((always_inline)) INLINE static void rt_init_part(
* @param cosmo Cosmology.
*/
__attribute__((always_inline)) INLINE static void rt_reset_part(
- struct part* restrict p, const struct cosmology* cosmo) {
+ struct part* restrict p, const struct cosmology* cosmo) {}
+
+/**
+ * @brief Reset RT particle data which needs to be reset each sub-cycle.
+ *
+ * @param p the particle to work on
+ * @param cosmo Cosmology.
+ * @param dt the current particle RT time step
+ */
+__attribute__((always_inline)) INLINE static void rt_reset_part_each_subcycle(
+ struct part* restrict p, const struct cosmology* cosmo, double dt) {
struct rt_part_data* rpd = &p->rt_data;
@@ -93,15 +103,7 @@ __attribute__((always_inline)) INLINE static void rt_reset_part(
rt_check_unphysical_state(&rpd->conserved[g].urad, rpd->conserved[g].frad,
urad_old, cred);
}
-}
-
-/**
- * @brief Reset RT particle data which needs to be reset each sub-cycle.
- *
- * @param p the particle to work on
- */
-__attribute__((always_inline)) INLINE static void rt_reset_part_each_subcycle(
- struct part* restrict p){};
+};
/**
* @brief First initialisation of the RT hydro particle data.
@@ -132,6 +134,7 @@ __attribute__((always_inline)) INLINE static void rt_first_init_part(
rt_init_part(p);
rt_reset_part(p, cosmo);
+ rt_reset_part_each_subcycle(p, cosmo, 0.);
}
/**
diff --git a/src/rt/debug/rt.h b/src/rt/debug/rt.h
index b36e0b2972d50608ec90a7db1e553417942de8a8..4b65951a634b830b9e3f23bfbf7a0dd4e7e52995 100644
--- a/src/rt/debug/rt.h
+++ b/src/rt/debug/rt.h
@@ -87,9 +87,11 @@ __attribute__((always_inline)) INLINE static void rt_reset_part(
* @brief Reset RT particle data which needs to be reset each sub-cycle.
*
* @param p the particle to work on
+ * @param cosmo Cosmology.
+ * @param dt the current particle RT time step
*/
__attribute__((always_inline)) INLINE static void rt_reset_part_each_subcycle(
- struct part* restrict p) {
+ struct part* restrict p, const struct cosmology* cosmo, double dt) {
rt_debugging_reset_each_subcycle(p);
}
@@ -107,7 +109,7 @@ __attribute__((always_inline)) INLINE static void rt_first_init_part(
rt_init_part(p);
rt_reset_part(p, cosmo);
- rt_reset_part_each_subcycle(p);
+ rt_reset_part_each_subcycle(p, cosmo, 0.);
p->rt_data.debug_radiation_absorbed_tot = 0ULL;
}
diff --git a/src/rt/none/rt.h b/src/rt/none/rt.h
index 5f256cfc83deeaf7d2980de79f44316a763d42d2..3555a4953eb90c1d7aff4c5505058693b7179ab6 100644
--- a/src/rt/none/rt.h
+++ b/src/rt/none/rt.h
@@ -76,9 +76,11 @@ __attribute__((always_inline)) INLINE static void rt_reset_part(
* @brief Reset RT particle data which needs to be reset each sub-cycle.
*
* @param p the particle to work on
+ * @param cosmo Cosmology.
+ * @param dt the current particle RT time step
*/
__attribute__((always_inline)) INLINE static void rt_reset_part_each_subcycle(
- struct part* restrict p) {}
+ struct part* restrict p, const struct cosmology* cosmo, double dt) {}
/**
* @brief First initialisation of the RT hydro particle data.
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index e44ad9a05e7c3b461a4662d9d01fa4e8ddd89907..ce87f19dbcc856cc7f9b28eb9d8ed660b0f992d5 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -1603,6 +1603,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
void runner_do_rt_ghost1(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology *cosmo = e->cosmology;
int count = c->hydro.count;
/* Anything to do here? */
@@ -1631,7 +1633,17 @@ void runner_do_rt_ghost1(struct runner *r, struct cell *c, int timer) {
/* First reset everything that needs to be reset for the following
* subcycle */
- rt_reset_part_each_subcycle(p);
+ const integertime_t ti_current_subcycle = e->ti_current_subcycle;
+ const integertime_t ti_step =
+ get_integer_timestep(p->rt_time_data.time_bin);
+ const integertime_t ti_begin = get_integer_time_begin(
+ ti_current_subcycle + 1, p->rt_time_data.time_bin);
+ const integertime_t ti_end = ti_begin + ti_step;
+
+ const float dt =
+ rt_part_dt(ti_begin, ti_end, e->time_base, with_cosmology, cosmo);
+
+ rt_reset_part_each_subcycle(p, cosmo, dt);
/* Now finish up injection */
rt_finalise_injection(p, e->rt_props);
diff --git a/src/runner_main.c b/src/runner_main.c
index 4ef7ddbfdbaf664872dea94d7e1d63b3061cb0e9..17a4bf1246a2553ef95953b02ffa19f3bcdc6801 100644
--- a/src/runner_main.c
+++ b/src/runner_main.c
@@ -520,7 +520,7 @@ void *runner_main(void *data) {
} else if (t->subtype == task_subtype_rt_gradient) {
runner_do_recv_part(r, ci, 2, 1);
} else if (t->subtype == task_subtype_rt_transport) {
- runner_do_recv_part(r, ci, 0, 1);
+ runner_do_recv_part(r, ci, -1, 1);
} else if (t->subtype == task_subtype_part_swallow) {
cell_unpack_part_swallow(ci,
(struct black_holes_part_data *)t->buff);
diff --git a/src/runner_recv.c b/src/runner_recv.c
index f225e1e6d6c7567b9cc06fbfccdae70e6b4e0017..c36f6d2349c29a6aa042f488ec9d49b524168153 100644
--- a/src/runner_recv.c
+++ b/src/runner_recv.c
@@ -70,6 +70,18 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int clear_sorts,
* need to force a sort after the gradient recv. */
clear_sorts = !cell_get_flag(c, cell_flag_skip_rt_sort);
cell_clear_flag(c, cell_flag_skip_rt_sort);
+ } else if (clear_sorts == -1) {
+ /* When running with RT and symmetric MPI exchanges, the first RT
+ * recv tasks (gradients) may not necessarily run, i.e. in cases
+ * where a local cell is inactive and needs to be updated by an
+ * active foreign cell. Should that occur, we need to make sure
+ * that the skip_rt_sort cell flag has been properly cleaned up
+ * first. */
+ cell_clear_flag(c, cell_flag_skip_rt_sort);
+ /* clear_sorts == -1 is abused to mark the case where this
+ * function is called for a recv_rt_transport task. We're not
+ * actually sorting now, so don't clear the sorts. */
+ clear_sorts = 0;
}
/* Clear this cell's sorted mask. */
diff --git a/src/scheduler.c b/src/scheduler.c
index 0e9c6f66132b7252134f0bb5b1b14593ca34574a..4a9a771c6c9466bc2c5ce02d90db1cfa713df250 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1968,8 +1968,7 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
t->weight = 0.f;
for (int j = 0; j < t->nr_unlock_tasks; j++)
- if (t->unlock_tasks[j]->weight > t->weight)
- t->weight = t->unlock_tasks[j]->weight;
+ t->weight += t->unlock_tasks[j]->weight;
const float count_i = (t->ci != NULL) ? t->ci->hydro.count : 0.f;
const float count_j = (t->cj != NULL) ? t->cj->hydro.count : 0.f;
@@ -2447,8 +2446,6 @@ void scheduler_start(struct scheduler *s) {
* @param t The #task.
*/
void scheduler_enqueue(struct scheduler *s, struct task *t) {
- /* The target queue for this task. */
- int qid = -1;
/* Ignore skipped tasks */
if (t->skip) return;
@@ -2482,16 +2479,21 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
#endif
/* Find the previous owner for each task type, and do
- any pre-processing needed. */
+ * any pre-processing needed. */
+ short int qid = -1;
+ short int *owner = NULL;
switch (t->type) {
case task_type_self:
case task_type_sub_self:
if (t->subtype == task_subtype_grav ||
t->subtype == task_subtype_grav_bkg ||
- t->subtype == task_subtype_external_grav)
+ t->subtype == task_subtype_external_grav) {
qid = t->ci->grav.super->owner;
- else
+ owner = &t->ci->grav.super->owner;
+ } else {
qid = t->ci->hydro.super->owner;
+ owner = &t->ci->hydro.super->owner;
+ }
break;
case task_type_sort:
case task_type_stars_resort:
@@ -2507,6 +2509,7 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
case task_type_drift_sink:
case task_type_cooling:
qid = t->ci->hydro.super->owner;
+ owner = &t->ci->hydro.super->owner;
break;
case task_type_grav_down:
case task_type_grav_long_range:
@@ -2519,6 +2522,7 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
case task_type_pack:
case task_type_unpack:
qid = t->ci->grav.super->owner;
+ owner = &t->ci->grav.super->owner;
break;
case task_type_kick1:
case task_type_kick2:
@@ -2527,25 +2531,16 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
case task_type_timestep_limiter:
case task_type_timestep_sync:
qid = t->ci->super->owner;
+ owner = &t->ci->super->owner;
break;
case task_type_pair:
case task_type_sub_pair:
- if (t->subtype == task_subtype_grav ||
- t->subtype == task_subtype_grav_zoombkg ||
- t->subtype == task_subtype_grav_bkg ||
- t->subtype == task_subtype_grav_bkgzoom ||
- t->subtype == task_subtype_external_grav) {
- qid = t->ci->grav.super->owner;
- if (qid < 0 ||
- s->queues[qid].count > s->queues[t->cj->grav.super->owner].count)
- qid = t->cj->grav.super->owner;
- } else if (t->subtype == task_subtype_grav_bkg_pool) {
- qid = t->ci->grav.super->owner;
- } else {
- qid = t->ci->hydro.super->owner;
- if (qid < 0 ||
- s->queues[qid].count > s->queues[t->cj->hydro.super->owner].count)
- qid = t->cj->hydro.super->owner;
+ qid = t->ci->super->owner;
+ owner = &t->ci->super->owner;
+ if (qid < 0 ||
+ s->queues[qid].count > s->queues[t->cj->super->owner].count) {
+ qid = t->cj->super->owner;
+ owner = &t->cj->super->owner;
}
break;
case task_type_recv:
@@ -2767,10 +2762,12 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
if (qid >= s->nr_queues) error("Bad computed qid.");
- /* If no previous owner, pick a random queue. */
- /* Note that getticks() is random enough */
- if (qid < 0) qid = getticks() % s->nr_queues;
-
+ /* If no qid, pick a random queue. */
+ if (qid < 0) qid = rand() % s->nr_queues;
+
+ /* Save qid as owner for next time a task accesses this cell. */
+ if (owner != NULL) *owner = qid;
+
/* Increase the waiting counter. */
atomic_inc(&s->waiting);
diff --git a/src/space.c b/src/space.c
index 970fd3338425350e686839439644b2b055ea6abc..55a9d06dc898d49b44f3a8aaa8152f98a2d5e4ac 100644
--- a/src/space.c
+++ b/src/space.c
@@ -474,59 +474,55 @@ void space_map_cells_pre(struct space *s, int full,
* @param nr_cells Number of #cell to pick up.
* @param cells Array of @c nr_cells #cell pointers in which to store the
* new cells.
+ * @param tpid ID of threadpool threadpool associated with cells_sub.
*/
void space_getcells(struct space *s, int nr_cells, struct cell **cells,
- const int thread_id) {
+ const short int tpid) {
/* For each requested cell... */
for (int j = 0; j < nr_cells; j++) {
/* Is the cell buffer empty? */
- if (s->cells_sub[thread_id] == NULL) {
-
- if (swift_memalign("cells_sub", (void **)&s->cells_sub[thread_id],
- cell_align,
+ if (s->cells_sub[tpid] == NULL) {
+ if (swift_memalign("cells_sub", (void **)&s->cells_sub[tpid], cell_align,
space_cellallocchunk * sizeof(struct cell)) != 0)
error("Failed to allocate more cells.");
/* Clear the newly-allocated cells. */
- bzero(s->cells_sub[thread_id],
- sizeof(struct cell) * space_cellallocchunk);
+ bzero(s->cells_sub[tpid], sizeof(struct cell) * space_cellallocchunk);
/* Constructed a linked list */
for (int k = 0; k < space_cellallocchunk - 1; k++)
- s->cells_sub[thread_id][k].next = &s->cells_sub[thread_id][k + 1];
- /* Signal the end of the list */
- s->cells_sub[thread_id][space_cellallocchunk - 1].next = NULL;
+ s->cells_sub[tpid][k].next = &s->cells_sub[tpid][k + 1];
+ s->cells_sub[tpid][space_cellallocchunk - 1].next = NULL;
}
/* Is the multipole buffer empty? */
- if (s->with_self_gravity && s->multipoles_sub[thread_id] == NULL) {
+ if (s->with_self_gravity && s->multipoles_sub[tpid] == NULL) {
if (swift_memalign(
- "multipoles_sub", (void **)&s->multipoles_sub[thread_id],
+ "multipoles_sub", (void **)&s->multipoles_sub[tpid],
multipole_align,
space_cellallocchunk * sizeof(struct gravity_tensors)) != 0)
error("Failed to allocate more multipoles.");
/* Constructed a linked list */
for (int k = 0; k < space_cellallocchunk - 1; k++)
- s->multipoles_sub[thread_id][k].next =
- &s->multipoles_sub[thread_id][k + 1];
- /* Signal the end of the list */
- s->multipoles_sub[thread_id][space_cellallocchunk - 1].next = NULL;
+ s->multipoles_sub[tpid][k].next = &s->multipoles_sub[tpid][k + 1];
+ s->multipoles_sub[tpid][space_cellallocchunk - 1].next = NULL;
}
/* Pick off the next cell. */
- cells[j] = s->cells_sub[thread_id];
- s->cells_sub[thread_id] = cells[j]->next;
+ cells[j] = s->cells_sub[tpid];
+ s->cells_sub[tpid] = cells[j]->next;
/* Hook the multipole */
if (s->with_self_gravity) {
- cells[j]->grav.multipole = s->multipoles_sub[thread_id];
- s->multipoles_sub[thread_id] = cells[j]->grav.multipole->next;
+ cells[j]->grav.multipole = s->multipoles_sub[tpid];
+ s->multipoles_sub[tpid] = cells[j]->grav.multipole->next;
}
}
+ /* Unlock the space. */
atomic_add(&s->tot_cells, nr_cells);
/* Init some things in the cell we just got. */
@@ -538,7 +534,7 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells,
bzero(cells[j], sizeof(struct cell));
cells[j]->grav.multipole = temp;
cells[j]->nodeID = -1;
- cells[j]->owner = thread_id;
+ cells[j]->tpid = tpid;
if (lock_init(&cells[j]->hydro.lock) != 0 ||
lock_init(&cells[j]->grav.plock) != 0 ||
lock_init(&cells[j]->grav.mlock) != 0 ||
@@ -558,9 +554,8 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells,
* @param s The #space.
*/
void space_free_buff_sort_indices(struct space *s) {
-
- for (int tid = 0; tid < s->e->nr_pool_threads; ++tid) {
- for (struct cell *finger = s->cells_sub[tid]; finger != NULL;
+ for (short int tpid = 0; tpid < s->e->nr_pool_threads; ++tpid) {
+ for (struct cell *finger = s->cells_sub[tpid]; finger != NULL;
finger = finger->next) {
cell_free_hydro_sorts(finger);
cell_free_stars_sorts(finger);
diff --git a/src/space.h b/src/space.h
index cc9253bb4b30bb74848cfa757656817cfcc516e2..f081521574ce6e7082168340391da44afd34c52a 100644
--- a/src/space.h
+++ b/src/space.h
@@ -159,13 +159,13 @@ struct space {
/*! The (level 0) cells themselves. */
struct cell *cells_top;
- /*! Buffer of unused cells for the sub-cells. */
+ /*! Buffer of unused cells for the sub-cells. One chunk per thread. */
struct cell **cells_sub;
/*! The multipoles associated with the top-level (level 0) cells */
struct gravity_tensors *multipoles_top;
- /*! Buffer of unused multipoles for the sub-cells. */
+ /*! Buffer of unused multipoles for the sub-cells. One chunk per thread. */
struct gravity_tensors **multipoles_sub;
/*! The indices of the *local* top-level cells */
@@ -499,7 +499,7 @@ void space_bparts_sort(struct bpart *bparts, int *ind, int *counts,
void space_sinks_sort(struct sink *sinks, int *ind, int *counts, int num_bins,
ptrdiff_t sinks_offset);
void space_getcells(struct space *s, int nr_cells, struct cell **cells,
- const int thread_id);
+ const short int tid);
void space_init(struct space *s, struct swift_params *params,
const struct cosmology *cosmo, double dim[3],
const struct hydro_props *hydro_properties,
diff --git a/src/space_recycle.c b/src/space_recycle.c
index a9f450bd7262d4853722b580e394dbf7bb5ae011..afa9cefbc3a4be93cbe30f6027a97eb2339710d8 100644
--- a/src/space_recycle.c
+++ b/src/space_recycle.c
@@ -255,28 +255,16 @@ void space_recycle(struct space *s, struct cell *c, const int lock) {
lock_destroy(&c->stars.star_formation_lock) != 0)
error("Failed to destroy spinlocks.");
- /* Lock the space. */
- if (lock) lock_lock(&s->lock);
-
- /* Thread which allocated this cell */
- const int owner = c->owner;
-
/* Hook the multipole back in the buffer */
if (s->with_self_gravity) {
- c->grav.multipole->next = s->multipoles_sub[owner];
- s->multipoles_sub[owner] = c->grav.multipole;
+ c->grav.multipole->next = s->multipoles_sub[c->tpid];
+ s->multipoles_sub[c->tpid] = c->grav.multipole;
}
/* Hook this cell into the buffer. */
- c->next = s->cells_sub[owner];
- s->cells_sub[owner] = c;
- if (lock)
- s->tot_cells -= 1;
- else
- atomic_dec(&s->tot_cells);
-
- /* Unlock the space. */
- if (lock) lock_unlock_blind(&s->lock);
+ c->next = s->cells_sub[c->tpid];
+ s->cells_sub[c->tpid] = c;
+ atomic_dec(&s->tot_cells);
}
/**
@@ -324,18 +312,17 @@ void space_recycle_list(struct space *s, struct cell *cell_list_begin,
/* Lock the space. */
lock_lock(&s->lock);
- /* Thread which allocated this cell */
- const int owner = cell_list_begin->owner;
-
- /* Hook the cells into the buffer. */
- cell_list_end->next = s->cells_sub[owner];
- s->cells_sub[owner] = cell_list_begin;
+ /* Hook the cells into the buffer keeping tpid if we can. */
+ short int tpid = cell_list_begin->tpid;
+ if (tpid < 0) tpid = 0;
+ cell_list_end->next = s->cells_sub[tpid];
+ s->cells_sub[tpid] = cell_list_begin;
atomic_sub(&s->tot_cells, count);
/* Hook the multipoles into the buffer. */
if (s->with_self_gravity) {
- multipole_list_end->next = s->multipoles_sub[owner];
- s->multipoles_sub[owner] = multipole_list_begin;
+ multipole_list_end->next = s->multipoles_sub[tpid];
+ s->multipoles_sub[tpid] = multipole_list_begin;
}
/* Unlock the space. */
diff --git a/src/space_split.c b/src/space_split.c
index e018cabbe29285ed2d654630f6b0b377c68677fd..811dff01dfc90a4cd50ce933db9396351c5944fe 100644
--- a/src/space_split.c
+++ b/src/space_split.c
@@ -57,7 +57,7 @@ void space_split_recursive(struct space *s, struct cell *c,
struct cell_buff *restrict bbuff,
struct cell_buff *restrict gbuff,
struct cell_buff *restrict sink_buff,
- const int thread_id) {
+ const short int tpid) {
const int count = c->hydro.count;
const int gcount = c->grav.count;
@@ -96,11 +96,9 @@ void space_split_recursive(struct space *s, struct cell *c,
struct engine *e = s->e;
const integertime_t ti_current = e->ti_current;
- /* Set the top level cell owner. Doing it here ensures top level cells
- * have the same owner as their progeny. */
- if (depth == 0) {
- c->owner = thread_id;
- }
+ /* Set the top level cell tpid. Doing it here ensures top level cells
+ * have the same tpid as their progeny. */
+ if (depth == 0) c->tpid = tpid;
/* If the buff is NULL, allocate it, and remember to free it. */
const int allocate_buffer = (buff == NULL && gbuff == NULL && sbuff == NULL &&
@@ -221,7 +219,7 @@ void space_split_recursive(struct space *s, struct cell *c,
c->split = 1;
/* Create the cell's progeny. */
- space_getcells(s, 8, c->progeny, thread_id);
+ space_getcells(s, 8, c->progeny, tpid);
for (int k = 0; k < 8; k++) {
struct cell *cp = c->progeny[k];
cp->hydro.count = 0;
@@ -719,6 +717,9 @@ void space_split_recursive(struct space *s, struct cell *c,
c->black_holes.h_max_active = black_holes_h_max_active;
c->maxdepth = maxdepth;
+ /* No runner owns this cell yet. We assign those during scheduling. */
+ c->owner = -1;
+
/* Store the global max depth */
if (c->depth == 0) atomic_max(&s->maxdepth, maxdepth);
@@ -773,10 +774,13 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data,
float max_softening = 0.f;
float max_mpole_power[SELF_GRAVITY_MULTIPOLE_ORDER + 1] = {0.f};
+ /* Threadpool id of current thread. */
+ short int tpid = threadpool_gettid();
+
/* Loop over the non-empty cells */
for (int ind = 0; ind < num_cells; ind++) {
struct cell *c = &cells_top[local_cells_with_particles[ind]];
- space_split_recursive(s, c, NULL, NULL, NULL, NULL, NULL, tid);
+ space_split_recursive(s, c, NULL, NULL, NULL, NULL, NULL, tpid);
if (s->with_self_gravity) {
min_a_grav =
diff --git a/src/threadpool.c b/src/threadpool.c
index c50b088e9c6077f579cea06307112616888c30c9..d1018fb89edc9d9c3d5f69d910d93e0e9ff3cc3b 100644
--- a/src/threadpool.c
+++ b/src/threadpool.c
@@ -39,6 +39,9 @@
#include "error.h"
#include "minmax.h"
+/* Keys for thread specific data. */
+static pthread_key_t threadpool_tid;
+
#ifdef SWIFT_DEBUG_THREADPOOL
/**
* @brief Store a log entry of the given chunk.
@@ -137,6 +140,10 @@ void threadpool_dump_log(struct threadpool *tp, const char *filename,
*/
static void threadpool_chomp(struct threadpool *tp, const int tid) {
+ /* Store the thread ID as thread specific data. */
+ int localtid = tid;
+ pthread_setspecific(threadpool_tid, &localtid);
+
/* Loop until we can't get a chunk. */
while (1) {
/* Compute the desired chunk size. */
@@ -164,14 +171,10 @@ static void threadpool_chomp(struct threadpool *tp, const int tid) {
#ifdef SWIFT_DEBUG_THREADPOOL
const ticks tic = getticks();
#endif
- if (tp->pass_tid) {
- tp->map_function_tid(
- (char *)tp->map_data + (tp->map_data_stride * task_ind), chunk_size,
- tp->map_extra_data, tid);
- } else {
- tp->map_function((char *)tp->map_data + (tp->map_data_stride * task_ind),
- chunk_size, tp->map_extra_data);
- }
+
+ tp->map_function((char *)tp->map_data + (tp->map_data_stride * task_ind),
+ chunk_size, tp->map_extra_data);
+
#ifdef SWIFT_DEBUG_THREADPOOL
threadpool_log(tp, tid, chunk_size, tic, getticks());
#endif
@@ -212,6 +215,13 @@ void threadpool_init(struct threadpool *tp, int num_threads) {
/* Initialize the thread counters. */
tp->num_threads = num_threads;
+ /* Create thread local data areas. Only do this once for all threads. */
+ pthread_key_create(&threadpool_tid, NULL);
+
+ /* Store the main thread ID as thread specific data. */
+ static int localtid = 0;
+ pthread_setspecific(threadpool_tid, &localtid);
+
#ifdef SWIFT_DEBUG_THREADPOOL
if ((tp->logs = (struct mapper_log *)malloc(sizeof(struct mapper_log) *
num_threads)) == NULL)
@@ -455,3 +465,11 @@ void threadpool_clean(struct threadpool *tp) {
free(tp->logs);
#endif
}
+
+/**
+ * @brief return the threadpool id of the current thread.
+ */
+int threadpool_gettid() {
+ int *tid = (int *)pthread_getspecific(threadpool_tid);
+ return *tid;
+}
diff --git a/src/threadpool.h b/src/threadpool.h
index 84d24a02eeaaded95e2cca991e7970fc1c31a23d..73e263157e65d9dbf3410142e676207a7a826c62 100644
--- a/src/threadpool.h
+++ b/src/threadpool.h
@@ -108,10 +108,7 @@ void threadpool_init(struct threadpool *tp, int num_threads);
void threadpool_map(struct threadpool *tp, threadpool_map_function map_function,
void *map_data, size_t N, int stride, int chunk,
void *extra_data);
-void threadpool_map_with_tid(struct threadpool *tp,
- threadpool_map_function_tid map_function,
- void *map_data, size_t N, int stride, int chunk,
- void *extra_data);
+int threadpool_gettid(void);
void threadpool_clean(struct threadpool *tp);
#ifdef SWIFT_DEBUG_THREADPOOL
void threadpool_reset_log(struct threadpool *tp);
diff --git a/theory/Gizmo/gizmo-implementation-details/bert.tex b/theory/Gizmo/gizmo-implementation-details/bert.tex
index 99975349d90ef99f736c37dd3a4ec8d2092189cb..455d7236f8ea2a18c3f8c22945240bec8a7f3256 100644
--- a/theory/Gizmo/gizmo-implementation-details/bert.tex
+++ b/theory/Gizmo/gizmo-implementation-details/bert.tex
@@ -78,7 +78,7 @@ The time derivatives are given by (eqn. (A4) in Hopkins)
\vec{w'} & \rho' & 0\\
0 & \vec{w'} & 1/\rho'\\
0 & \gamma P' & \vec{w'}
- \end{pmatrix} \vec{\nabla}X_k,
+ \end{pmatrix} \vec{\nabla}X_k, \label{eq:euler_primitive}
\end{equation}
or, for example~:
\begin{equation}
@@ -101,4 +101,47 @@ with $\vec{n}_k$ the unit vector along the coordinate axes ($k=x,y,z$) and $e =
Finally, we use the fluxes to update the conserved variables (eqn. (23) in Hopkins)~:
\begin{equation}
\Delta C_k = -\Delta t \sum_j \vec{F_{ij}}.\vec{A_{ij}}.
-\end{equation}
\ No newline at end of file
+\end{equation}
+
+\subsection{Other equations that are used}
+
+In order to couple the GIZMO schemes to subgrid physics, we need to provide a way to obtain the time
+derivative of the internal energy, which is not a basic quantity in these schemes. This time derivative can be
+obtained from the Euler equations, using the gradients. From \eqref{eq:euler_primitive}, we get
+
+\begin{equation}
+\frac{\partial{}P}{\partial{}t} = -\gamma{}P \vec{\nabla{}}.\vec{w} - \vec{w}.\vec{\nabla{}}P,
+\end{equation}
+while we also have
+\begin{equation}
+\frac{\partial{}P}{\partial{}t} = \frac{\partial{}}{\partial{}t} \left(
+ (\gamma{}-1) \rho{} u
+\right) = (\gamma{}-1)\rho{}\frac{\partial{}u}{\partial{}t} + (\gamma{}-1)u\frac{\partial{}\rho{}}{\partial{}t}
+\end{equation}
+or hence
+\begin{equation}
+\frac{\partial{}u}{\partial{}t} = \frac{1}{(\gamma{}-1)\rho{}} \frac{\partial{}P}{\partial{}t} -
+\frac{u}{\rho{}} \left( -\vec{w}.\vec{\nabla{}}\rho{} - \rho{}\vec{\nabla{}}.\vec{w} \right),
+\end{equation}
+where we again used \eqref{eq:euler_primitive} in the last term.
+
+Combining these equations, we get
+\begin{align}
+\frac{\partial{}u}{\partial{}t} &= \frac{1}{(\gamma{}-1)\rho{}} \left(
+-\gamma{}P \vec{\nabla{}}.\vec{w} - \vec{w}.\vec{\nabla{}}P
+\right) -
+\frac{u}{\rho{}} \left( -\vec{w}.\vec{\nabla{}}\rho{} - \rho{}\vec{\nabla{}}.\vec{w} \right) \\
+&= -\gamma{}u \vec{\nabla{}}.\vec{w} - \frac{1}{(\gamma{}-1)\rho{}} \vec{w}.\vec{\nabla{}}P
++ \frac{u}{\rho{}}\vec{w}.\vec{\nabla{}}\rho{} + u \vec{\nabla{}}.\vec{w} \\
+&= -(\gamma{}-1)u \vec{\nabla{}}.\vec{w} - \frac{1}{(\gamma{}-1)\rho{}} \vec{w}.\vec{\nabla{}}P
++ \frac{u}{\rho{}}\vec{w}.\vec{\nabla{}}\rho{}.
+\end{align}
+
+For convenience, we can eliminate $u$ from this equation again:
+\begin{equation}
+\frac{\partial{}u}{\partial{}t} = -\frac{P}{\rho{}}\vec{\nabla{}}.\vec{w}
+- \frac{1}{(\gamma{}-1)\rho{}} \vec{w} . \left(
+\vec{\nabla{}}P - \frac{P}{\rho{}} \vec{\nabla{}}\rho{}
+\right),
+\end{equation}
+where we recognise a factor $\vec{\nabla{}}u$ in the second term.
diff --git a/tools/task_plots/analyse_threadpool_tasks.py b/tools/task_plots/analyse_threadpool_tasks.py
index 5c25b7673bde02bb950f332d74382e2e8dc2e85b..34f75ea0732470732953fe40b0ce95ef7512c7ef 100755
--- a/tools/task_plots/analyse_threadpool_tasks.py
+++ b/tools/task_plots/analyse_threadpool_tasks.py
@@ -28,11 +28,6 @@ You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
-import matplotlib
-
-matplotlib.use("Agg")
-import matplotlib.collections as collections
-import matplotlib.ticker as plticker
import pylab as pl
import sys
import argparse
diff --git a/tools/task_plots/iplot_tasks.py b/tools/task_plots/iplot_tasks.py
index a96ad977c2604581fa1c33a234dd1f19d4b0a62b..50e65ebf9b7b285ce6f8c55f20d1006db1c30e67 100755
--- a/tools/task_plots/iplot_tasks.py
+++ b/tools/task_plots/iplot_tasks.py
@@ -458,6 +458,35 @@ class Container:
else:
fig.canvas.mpl_connect("button_press_event", self.onclick)
+ # Space bar to dump all tasks. Use with caution...
+ fig.canvas.mpl_connect("key_press_event", self.dump)
+
+ def dump(self, event):
+ # Dump all tasks to the console sorted by tic.
+ xlow = float(event.inaxes.viewLim.x0)
+ xhigh = float(event.inaxes.viewLim.x1)
+
+ if event.key == " ":
+ dumps = {}
+ for thread in range(nthread):
+ tics = self.ltics[thread]
+ tocs = self.ltocs[thread]
+ labels = self.llabels[thread]
+ for i in range(len(tics)):
+ if (tics[i] > xlow and tics[i] < xhigh) or (
+ tocs[i] > xlow and tocs[i] < xhigh
+ ):
+ tic = "{0:.3f}".format(tics[i])
+ toc = "{0:.3f}".format(tocs[i])
+ dumps[tics[i]] = (
+ labels[i] + ", tic/toc = " + tic + " / " + toc
+ )
+ print("")
+ print("Tasks in time range: " + str(xlow) + " -> " + str(xhigh))
+ for key in sorted(dumps):
+ print(dumps[key])
+ print("")
+
def onclick(self, event):
# Find thread, then scan for bounded task.
try:
diff --git a/tools/task_plots/iplot_threadpool.py b/tools/task_plots/iplot_threadpool.py
new file mode 100755
index 0000000000000000000000000000000000000000..70ca679e0a22a87c45f710d47b4f915a4f3a50ce
--- /dev/null
+++ b/tools/task_plots/iplot_threadpool.py
@@ -0,0 +1,390 @@
+#!/usr/bin/env python3
+"""
+Interactive plot of a threadpool dump.
+
+Usage:
+ iplot_threadpool.py [options] input.dat
+
+where input.dat is a threadpool info file for a step. Use the '-Y interval'
+flag of the swift or swift_mpi commands to create these (these will need to be
+built with the --enable-threadpool-debugging configure option).
+
+The plot can be scrolled and zoomed using the standard matplotlib
+controls, the type of task at a point can be queried by a mouse click
+(unless the --motion option is in effect when a continuous readout is
+shown) the task type and tic/toc range are reported in the terminal.
+
+Requires the tkinter module.
+
+This file is part of SWIFT.
+
+Copyright (C) 2022 Peter W. Draper (p.w.draper@durham.ac.uk)
+All Rights Reserved.
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import matplotlib
+
+matplotlib.use("TkAgg")
+import numpy as np
+import matplotlib.backends.backend_tkagg as tkagg
+from matplotlib.figure import Figure
+import tkinter as tk
+import math
+import matplotlib.collections as collections
+import matplotlib.ticker as plticker
+import pylab as pl
+import sys
+import argparse
+
+# Handle the command line.
+parser = argparse.ArgumentParser(description="Plot threadpool function graphs")
+
+parser.add_argument("input", help="Threadpool data file (-Y output)")
+parser.add_argument(
+ "-m",
+ "--motion",
+ dest="motion",
+ help="Track mouse motion, otherwise clicks (def: clicks)",
+ default=False,
+ action="store_true",
+)
+parser.add_argument(
+ "-l",
+ "--limit",
+ dest="limit",
+ help="Upper time limit in millisecs (def: depends on data)",
+ default=0,
+ type=float,
+)
+parser.add_argument(
+ "--height",
+ dest="height",
+ help="Height of plot in inches (def: 4)",
+ default=4.0,
+ type=float,
+)
+parser.add_argument(
+ "--width",
+ dest="width",
+ help="Width of plot in inches (def: 16)",
+ default=16.0,
+ type=float,
+)
+parser.add_argument(
+ "-v",
+ "--verbose",
+ dest="verbose",
+ help="Show colour assignments and other details (def: False)",
+ default=False,
+ action="store_true",
+)
+
+args = parser.parse_args()
+infile = args.input
+delta_t = args.limit
+
+# Basic plot configuration.
+PLOT_PARAMS = {
+ "axes.labelsize": 10,
+ "axes.titlesize": 10,
+ "font.size": 12,
+ "xtick.labelsize": 10,
+ "ytick.labelsize": 10,
+ "figure.figsize": (args.width, args.height),
+ "figure.subplot.left": 0.03,
+ "figure.subplot.right": 0.995,
+ "figure.subplot.bottom": 0.09,
+ "figure.subplot.top": 0.99,
+ "figure.subplot.wspace": 0.0,
+ "figure.subplot.hspace": 0.0,
+ "lines.markersize": 6,
+ "lines.linewidth": 3.0,
+}
+pl.rcParams.update(PLOT_PARAMS)
+
+# A number of colours for the various types. Recycled when there are
+# more task types than colours...
+colours = [
+ "cyan",
+ "lightgray",
+ "darkblue",
+ "yellow",
+ "tan",
+ "dodgerblue",
+ "sienna",
+ "aquamarine",
+ "bisque",
+ "blue",
+ "green",
+ "lightgreen",
+ "brown",
+ "purple",
+ "moccasin",
+ "olivedrab",
+ "chartreuse",
+ "olive",
+ "darkgreen",
+ "green",
+ "mediumseagreen",
+ "mediumaquamarine",
+ "darkslategrey",
+ "mediumturquoise",
+ "black",
+ "cadetblue",
+ "skyblue",
+ "red",
+ "slategray",
+ "gold",
+ "slateblue",
+ "blueviolet",
+ "mediumorchid",
+ "firebrick",
+ "magenta",
+ "hotpink",
+ "pink",
+ "orange",
+ "lightgreen",
+]
+maxcolours = len(colours)
+
+# Read header. First two lines.
+with open(infile) as infid:
+ head = [next(infid) for x in range(2)]
+header = head[1][2:].strip()
+header = eval(header)
+nthread = int(header["num_threads"]) + 1
+CPU_CLOCK = float(header["cpufreq"]) / 1000.0
+print("Number of threads: ", nthread)
+if args.verbose:
+ print("CPU frequency:", CPU_CLOCK * 1000.0)
+
+# Read input.
+data = pl.genfromtxt(infile, dtype=None, delimiter=" ", encoding=None)
+
+# Mixed types, so need to separate.
+tics = []
+tocs = []
+funcs = []
+threads = []
+chunks = []
+for i in data:
+ if i[0] != "#":
+ funcs.append(i[0].replace("_mapper", ""))
+ if i[1] < 0:
+ threads.append(nthread - 1)
+ else:
+ threads.append(i[1])
+ chunks.append(i[2])
+ tics.append(i[3])
+ tocs.append(i[4])
+tics = pl.array(tics)
+tocs = pl.array(tocs)
+funcs = pl.array(funcs)
+threads = pl.array(threads)
+chunks = pl.array(chunks)
+
+# Recover the start and end time
+mintic_step = min(tics)
+tic_step = mintic_step
+toc_step = max(tocs)
+print("# Min tic = ", mintic_step)
+
+# Calculate the time range, if not given.
+delta_t = delta_t * CPU_CLOCK
+if delta_t == 0:
+ dt = toc_step - tic_step
+ if dt > delta_t:
+ delta_t = dt
+ print("Data range: ", delta_t / CPU_CLOCK, "ms")
+
+# Once more doing the real gather and plots this time.
+start_t = float(tic_step)
+tics -= tic_step
+tocs -= tic_step
+end_t = (toc_step - start_t) / CPU_CLOCK
+
+# Get all "task" names and assign colours.
+TASKTYPES = pl.unique(funcs)
+print(TASKTYPES)
+
+# Set colours of tasks.
+TASKCOLOURS = {}
+ncolours = 0
+for task in TASKTYPES:
+ TASKCOLOURS[task] = colours[ncolours]
+ ncolours = (ncolours + 1) % maxcolours
+
+# For fiddling with colours...
+if args.verbose:
+ print("#Selected colours:")
+ for task in sorted(TASKCOLOURS.keys()):
+ print("# " + task + ": " + TASKCOLOURS[task])
+
+tasks = {}
+tasks[-1] = []
+for i in range(nthread):
+ tasks[i] = []
+
+for i in range(len(threads)):
+ thread = threads[i]
+ tasks[thread].append({})
+ tasks[thread][-1]["type"] = funcs[i]
+ tic = tics[i] / CPU_CLOCK
+ toc = tocs[i] / CPU_CLOCK
+ tasks[thread][-1]["tic"] = tic
+ tasks[thread][-1]["toc"] = toc
+ tasks[thread][-1]["colour"] = TASKCOLOURS[funcs[i]]
+
+# Do the plotting.
+fig = Figure()
+ax = fig.add_subplot(1, 1, 1)
+ax.set_xlim(-delta_t * 0.01 / CPU_CLOCK, delta_t * 1.01 / CPU_CLOCK)
+ax.set_ylim(0.5, nthread + 1.0)
+
+ltics = []
+ltocs = []
+llabels = []
+for i in range(nthread):
+
+ # Collect ranges and colours into arrays. Also indexed lists for lookup tables.
+ tictocs = []
+ colours = []
+ tics = []
+ tocs = []
+ labels = []
+ for task in tasks[i]:
+ tictocs.append((task["tic"], task["toc"] - task["tic"]))
+ colours.append(task["colour"])
+
+ tics.append(task["tic"])
+ tocs.append(task["toc"])
+ labels.append(task["type"])
+
+ # Add to look up tables.
+ ltics.append(tics)
+ ltocs.append(tocs)
+ llabels.append(labels)
+
+ # Now plot.
+ ax.broken_barh(tictocs, [i + 0.55, 0.9], facecolors=colours, linewidth=0)
+
+# Start and end of time-step
+ax.plot([0, 0], [0, nthread + 1], "k--", linewidth=1)
+ax.plot([end_t, end_t], [0, nthread + 1], "k--", linewidth=1)
+
+# Labels.
+ax.set_xlabel("Wall clock time [ms]")
+ax.set_ylabel("Thread ID")
+
+loc = plticker.MultipleLocator(base=1)
+ax.yaxis.set_major_locator(loc)
+ax.grid(True, which="major", axis="y", linestyle="-")
+
+
+class Container:
+ def __init__(self, window, figure, motion, nthread, ltics, ltocs, llabels):
+ self.window = window
+ self.figure = figure
+ self.motion = motion
+ self.nthread = nthread
+ self.ltics = ltics
+ self.ltocs = ltocs
+ self.llabels = llabels
+
+ def plot(self):
+ canvas = tkagg.FigureCanvasTkAgg(self.figure, master=self.window)
+ wcanvas = canvas.get_tk_widget()
+ wcanvas.config(width=1000, height=300)
+ wcanvas.pack(side=tk.TOP, expand=True, fill=tk.BOTH)
+
+ toolbar = tkagg.NavigationToolbar2Tk(canvas, self.window)
+ toolbar.update()
+ self.output = tk.StringVar()
+ label = tk.Label(
+ self.window, textvariable=self.output, bg="white", fg="red", bd=2
+ )
+ label.pack(side=tk.RIGHT, expand=True, fill=tk.X)
+ wcanvas.pack(side=tk.TOP, expand=True, fill=tk.BOTH)
+
+ canvas.draw()
+
+ # Print task type using mouse clicks or motion.
+ if self.motion:
+ fig.canvas.mpl_connect("motion_notify_event", self.onclick)
+ else:
+ fig.canvas.mpl_connect("button_press_event", self.onclick)
+
+ # Space bar to dump all tasks. Use with caution...
+ fig.canvas.mpl_connect("key_press_event", self.dump)
+
+ def dump(self, event):
+ # Dump all tasks to the console sorted by tic.
+ xlow = float(event.inaxes.viewLim.x0)
+ xhigh = float(event.inaxes.viewLim.x1)
+
+ if event.key == " ":
+ dumps = {}
+ for thread in range(nthread):
+ tics = self.ltics[thread]
+ tocs = self.ltocs[thread]
+ labels = self.llabels[thread]
+ for i in range(len(tics)):
+ if (tics[i] > xlow and tics[i] < xhigh) or (
+ tocs[i] > xlow and tocs[i] < xhigh
+ ):
+ tic = "{0:.3f}".format(tics[i])
+ toc = "{0:.3f}".format(tocs[i])
+ dumps[tics[i]] = (
+ labels[i] + ", tic/toc = " + tic + " / " + toc
+ )
+ print("")
+ print("Tasks in time range: " + str(xlow) + " -> " + str(xhigh))
+ for key in sorted(dumps):
+ print(dumps[key])
+ print("")
+
+ def onclick(self, event):
+ # Find thread, then scan for bounded task.
+ try:
+ thread = int(round(event.ydata)) - 1
+ outstr = "none"
+ if thread >= 0 and thread < self.nthread:
+ tics = self.ltics[thread]
+ tocs = self.ltocs[thread]
+ labels = self.llabels[thread]
+ for i in range(len(tics)):
+ if event.xdata > tics[i] and event.xdata < tocs[i]:
+ tic = "{0:.3f}".format(tics[i])
+ toc = "{0:.3f}".format(tocs[i])
+ outstr = labels[i] + ", tic/toc = " + tic + " / " + toc
+ break
+ self.output.set(outstr)
+ except TypeError:
+ # Out of bounds clears field.
+ self.output.set("")
+ pass
+
+ def quit(self):
+ self.window.destroy()
+
+
+window = tk.Tk()
+window.protocol("WM_DELETE_WINDOW", window.quit)
+container = Container(window, fig, args.motion, nthread, ltics, ltocs, llabels)
+container.plot()
+window.mainloop()
+
+sys.exit(0)
diff --git a/tools/task_plots/plot_threadpool.py b/tools/task_plots/plot_threadpool.py
index 05f07a7939f1f28fb5b7d653e6e1b57db9776885..1da8fb2ac010c84041628c51fd7ff213aef06221 100755
--- a/tools/task_plots/plot_threadpool.py
+++ b/tools/task_plots/plot_threadpool.py
@@ -30,14 +30,15 @@ You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
+import argparse
+import sys
+import pylab as pl
+import matplotlib.ticker as plticker
+import matplotlib.collections as collections
+import math
import matplotlib
matplotlib.use("Agg")
-import matplotlib.collections as collections
-import matplotlib.ticker as plticker
-import pylab as pl
-import sys
-import argparse
# Handle the command line.
parser = argparse.ArgumentParser(description="Plot threadpool function graphs")
@@ -289,7 +290,8 @@ j = 0
for task in tasks[nthread - expand]:
tictocs.append((task["tic"], task["toc"] - task["tic"]))
colours.append(task["colour"])
-ax.broken_barh(tictocs, [0, (nthread - 1)], facecolors=colours, linewidth=0, alpha=0.15)
+ax.broken_barh(tictocs, [0, (nthread - 1)],
+ facecolors=colours, linewidth=0, alpha=0.15)
# And we don't plot the fake thread.
nthread = nthread - expand
@@ -313,19 +315,24 @@ for i in range(nthread):
ax.broken_barh(tictocs, [i + 0.05, 0.90], facecolors=colours, linewidth=0)
# Legend and room for it.
-nrow = len(typesseen) / 5
+nrow = math.ceil(len(typesseen) / 4)
if not args.nolegend:
- ax.fill_between([0, 0], nthread + 0.5, nthread + nrow + 0.5, facecolor="white")
+ ax.fill_between([0, 0], nthread, nthread + nrow, facecolor="white")
ax.set_ylim(0, nthread + 0.5)
ax.legend(
- loc="upper center", shadow=True, bbox_to_anchor=(0.0, 1.3, 1.0, 0.2), ncol=6, fontsize=9
+ loc="lower left",
+ shadow=True,
+ bbox_to_anchor=(0.0, 1.0, 1.0, 0.2),
+ mode="expand",
+ ncol=4,
)
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width, box.height * 0.65])
# Start and end of time-step
real_start_t = (mintic_step - tic_step) / CPU_CLOCK
-ax.plot([real_start_t, real_start_t], [0, nthread + nrow + 1], "k--", linewidth=1)
+ax.plot([real_start_t, real_start_t], [
+ 0, nthread + nrow + 1], "k--", linewidth=1)
ax.plot([end_t, end_t], [0, nthread + nrow + 1], "k--", linewidth=1)
@@ -340,8 +347,7 @@ loc = plticker.MultipleLocator(base=expand)
ax.yaxis.set_major_locator(loc)
ax.grid(True, which="major", axis="y", linestyle="-")
-pl.show()
-pl.savefig(outpng)
+pl.savefig(outpng, bbox_inches="tight", dpi=100)
print("Graphics done, output written to", outpng)
sys.exit(0)
diff --git a/tools/task_plots/process_plot_threadpool b/tools/task_plots/process_plot_threadpool
index b1238b3668067f17b98099f6287fa63fec5d2395..87b8726924a19122144df4aad3a1b5ac5cf4842c 100755
--- a/tools/task_plots/process_plot_threadpool
+++ b/tools/task_plots/process_plot_threadpool
@@ -59,12 +59,12 @@ fi
list=""
for f in $files; do
s=$(echo $f| sed 's,threadpool_info-step\(.*\).dat,\1,')
- list="$list $f $s poolstep${s}r"
+ list="$list $f $s threadpool-step${s}"
done
# And process them,
echo "Processing threadpool info files..."
-echo $list | xargs -P $NPROCS -n 3 /bin/bash -c "${SCRIPTHOME}/plot_threadpool.py --expand 1 --limit $TIMERANGE --width 16 --height 4 \$0 \$2 "
+echo $list | xargs -P $NPROCS -n 3 /bin/bash -c "${SCRIPTHOME}/plot_threadpool.py --expand 1 --limit $TIMERANGE --width 16 --height 8 \$0 \$2 "
echo $list | xargs -P $NPROCS -n 3 /bin/bash -c "${SCRIPTHOME}/analyse_threadpool_tasks.py \$0 --html > \$2.stats"
echo "Writing output threadpool-index.html file"
@@ -84,17 +84,19 @@ echo $list | xargs -n 3 | while read f s g; do
<h2>Step $s</h2>
EOF
cat <<EOF >> threadpool-index.html
-<a href="poolstep${s}r${i}.html"><img src="poolstep${s}r${i}.png" width=400px/></a>
+<a href="threadpool-step${s}.html"><img src="threadpool-step${s}.png" width=400px/></a>
EOF
- cat <<EOF > poolstep${s}r${i}.html
+ cat <<EOF > threadpool-step${s}.html
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<body>
-<img src="poolstep${s}r${i}.png">
+<img src="threadpool-step${s}.png">
<pre>
+<nav>Jump to: <a href="#all">all threads</a> <a href="#dead">dead times</a></nav>
EOF
-cat poolstep${s}r${i}.stats >> poolstep${s}r${i}.html
-cat <<EOF >> poolstep${s}r${i}.html
+cat threadpool-step${s}.stats >> threadpool-step${s}.html
+cat <<EOF >> threadpool-step${s}.html
+</pre>
</body>
</html>
EOF
diff --git a/tools/task_plots/process_plot_threadpool_MPI b/tools/task_plots/process_plot_threadpool_MPI
new file mode 100755
index 0000000000000000000000000000000000000000..826c60497ba17308aa7ff46cf82344475ba8894e
--- /dev/null
+++ b/tools/task_plots/process_plot_threadpool_MPI
@@ -0,0 +1,116 @@
+#!/bin/bash
+#
+# Usage:
+# process_plot_threadpool_MPI nprocess [time-range-ms]
+#
+# Description:
+# Process all the threadpool info files in the current directory
+# creating function graphs for steps and threads. MPI version
+#
+# The input files are created by a run using the "-Y interval" flag and
+# should be named "threadpool_info-rank<m>-step<n>.dat" in the current
+# directory. All located files will be processed using "nprocess" concurrent
+# processes and all plots will have the same time range if one is given.
+# An output HTML file "index.html" will be created to view all the plots.
+#
+#
+# This file is part of SWIFT:
+#
+# Copyright (C) 2022 Peter W. Draper (p.w.draper@durham.ac.uk)
+# All Rights Reserved.
+#
+# 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/>.
+
+# Handle command-line
+if test "$1" = ""; then
+ echo "Usage: $0 nprocess [time-range-ms]"
+ exit 1
+fi
+NPROCS=$1
+TIMERANGE=0
+LIMIT="(autoranged)"
+if test "$2" != ""; then
+ TIMERANGE=$2
+ LIMIT=""
+fi
+
+# Locate script.
+SCRIPTHOME=$(dirname "$0")
+
+# Find all thread info files. Use sort to get into correct order, step then rank.
+files=$(ls -1 threadpool_info-rank*-step*.dat| sort -t- -k3,3)
+if test $? != 0; then
+ echo "Failed to find any threadpool info files"
+ exit 1
+fi
+
+# Construct list of names, the step no and names for the graphics.
+list=""
+for f in $files; do
+ r=$(echo $f| sed 's,threadpool_info-rank\(.*\)-step.*.dat,\1,')
+ s=$(echo $f| sed 's,threadpool_info-rank.*-step\(.*\).dat,\1,')
+ list="$list $f $s $r threadpool-step${s}-rank${r}"
+done
+
+# And process them,
+echo "Processing threadpool info files..."
+echo $list | xargs -P $NPROCS -n 4 /bin/bash -c "${SCRIPTHOME}/plot_threadpool.py --expand 1 --limit $TIMERANGE --width 16 --height 4 \$0 \$3 "
+echo $list | xargs -P $NPROCS -n 4 /bin/bash -c "${SCRIPTHOME}/analyse_threadpool_tasks.py \$0 --html > \$3.stats"
+
+echo "Writing output threadpool-index.html file"
+# Construct document - serial.
+cat <<EOF > threadpool-index.html
+ <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+ <head>
+ <title>SWIFT threadpool tasks $LIMIT</title>
+ </head>
+ <body>
+ <h1>SWIFT threadpool tasks $LIMIT</h1>
+EOF
+
+echo $list | xargs -n 4 | while read f s r g; do
+ cat <<EOF >> threadpool-index.html
+<h2>Step $s Rank $r</h2>
+EOF
+ cat <<EOF >> threadpool-index.html
+<a href="threadpool-step${s}-rank${r}.html"><img src="threadpool-step${s}-rank${r}.png" width=400px/></a>
+EOF
+ cat <<EOF > threadpool-step${s}-rank${r}.html
+ <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<body>
+<img src="threadpool-step${s}-rank${r}.png">
+<pre>
+<nav>Jump to: <a href="#all">all threads</a> <a href="#dead">dead times</a></nav>
+
+<pre>
+EOF
+cat threadpool-step${s}-rank${r}.stats >> threadpool-step${s}-rank${r}.html
+cat <<EOF >> threadpool-step${s}-rank${r}.html
+</pre>
+</body>
+</html>
+EOF
+
+done
+
+cat <<EOF >> threadpool-index.html
+ </body>
+</html>
+EOF
+
+echo "Finished"
+
+exit