diff --git a/.gitignore b/.gitignore
index 1fd274009591c7ac98d71e05d536c98d4c17c485..9772f54a40e1663bf5441cb5fc8cddf03d22f127 100644
--- a/.gitignore
+++ b/.gitignore
@@ -30,6 +30,7 @@ examples/fof_mpi
examples/*/*/*.xmf
examples/*/*/*.dat
examples/*/*/*.png
+examples/*/*/*.pdf
examples/*/*/*.mp4
examples/*/*/*.txt
examples/*/*/*.rst
diff --git a/examples/Cooling/FeedbackEvent_3D/.gitignore b/examples/Cooling/FeedbackEvent_3D/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..3a8485fdcf1c2638f90686598a442bfef6f50f64
--- /dev/null
+++ b/examples/Cooling/FeedbackEvent_3D/.gitignore
@@ -0,0 +1,2 @@
+nocool_*/*
+default_*/*
diff --git a/examples/Cooling/FeedbackEvent_3D/README.md b/examples/Cooling/FeedbackEvent_3D/README.md
index 236d03484c0a0afe846b464aef64573788fc639e..565990dd302ccb334c67ffa234294c031b8e6d9f 100644
--- a/examples/Cooling/FeedbackEvent_3D/README.md
+++ b/examples/Cooling/FeedbackEvent_3D/README.md
@@ -16,4 +16,12 @@ This test emulates what the EAGLE model does to particles for feedback, i.e.
+ Heats a single particle to 10^7.5 K
+ Does _not_ switch off cooling
-+ Runs to completion.
\ No newline at end of file
++ Runs to completion.
+
+
+Running Multiple Tests
+----------------------
+
+If you would like to run a suite of tests, try the runs.sh script. You'll
+need to set the directories in the parameter file to be one higher, i.e.
+../coolingtables rather than ./coolingtables.
diff --git a/examples/Cooling/FeedbackEvent_3D/feedback.yml b/examples/Cooling/FeedbackEvent_3D/feedback.yml
index da624834dcd3347593f5d48c849a8c434cd70392..5e4ddfc618ca64e6633772f155abbdc021cedc0e 100644
--- a/examples/Cooling/FeedbackEvent_3D/feedback.yml
+++ b/examples/Cooling/FeedbackEvent_3D/feedback.yml
@@ -9,20 +9,20 @@ InternalUnitSystem:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1e-3 # The end time of the simulation (in internal units).
+ time_end: 1e-4 # The end time of the simulation (in internal units).
dt_min: 1e-9 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: feedback # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 1e-4 # Time difference between consecutive outputs (in internal units)
+ delta_time: 1e-5 # Time difference between consecutive outputs (in internal units)
compression: 1
# Parameters governing the conserved quantities statistics
Statistics:
- delta_time: 1e-3 # Time between statistics output
+ delta_time: 1e-6 # Time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
@@ -75,4 +75,16 @@ EAGLECooling:
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
\ No newline at end of file
+ He_reion_eV_p_H: 2.0
+
+# Parameters for the EAGLE "equation of state"
+EAGLEEntropyFloor:
+ Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+ Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+ Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+ Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+ Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+ Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
+
diff --git a/examples/Cooling/FeedbackEvent_3D/makeIC.py b/examples/Cooling/FeedbackEvent_3D/makeIC.py
index 0002dc459730c53500ab75e88d8765701e4937a2..24052cee995675742808df6cc16ae4b523389351 100644
--- a/examples/Cooling/FeedbackEvent_3D/makeIC.py
+++ b/examples/Cooling/FeedbackEvent_3D/makeIC.py
@@ -30,43 +30,44 @@ import numpy as np
# Parameters
gamma = 5.0 / 3.0
initial_density = 0.1 * mh / (cm ** 3)
-initial_temperature = 1e4 * K
+initial_temperature = 2550 * (5/4) * K # Equilibrium temperature at solar abundance
inject_temperature = 10 ** (7.5) * K
mu = 0.5
particle_mass = 1e6 * msun
-unit_system = cosmo_units
-file_name = "feedback.hdf5"
+if __name__ == "__main__":
+ unit_system = cosmo_units
+ file_name = "feedback.hdf5"
-# Read in glass file
-with h5py.File("glassCube_32.hdf5", "r") as handle:
- positions = handle["/PartType0/Coordinates"][:]
- h = handle["PartType0/SmoothingLength"][:] * 0.3
+ # Read in glass file
+ with h5py.File("glassCube_32.hdf5", "r") as handle:
+ positions = handle["/PartType0/Coordinates"][:]
+ h = handle["PartType0/SmoothingLength"][:] * 0.3
-number_of_particles = len(h)
-side_length = (number_of_particles * particle_mass / initial_density) ** (1 / 3)
-side_length.convert_to_base(unit_system)
+ number_of_particles = len(h)
+ side_length = (number_of_particles * particle_mass / initial_density) ** (1 / 3)
+ side_length.convert_to_base(unit_system)
-print(f"Your box has a side length of {side_length}")
+ print(f"Your box has a side length of {side_length}")
-# Find the central particle
-central_particle = np.sum((positions - 0.5) ** 2, axis=1).argmin()
+ # Find the central particle
+ central_particle = np.sum((positions - 0.5) ** 2, axis=1).argmin()
-# Inject the feedback into that central particle
-background_internal_energy = (
- (1.0 / (mu * mh)) * (kb / (gamma - 1.0)) * initial_temperature
-)
-heated_internal_energy = (1.0 / (mu * mh)) * (kb / (gamma - 1.0)) * inject_temperature
-internal_energy = np.ones_like(h) * background_internal_energy
-internal_energy[central_particle] = heated_internal_energy
+ # Inject the feedback into that central particle
+ background_internal_energy = (
+ (1.0 / (mu * mh)) * (kb / (gamma - 1.0)) * initial_temperature
+ )
+ heated_internal_energy = (1.0 / (mu * mh)) * (kb / (gamma - 1.0)) * inject_temperature
+ internal_energy = np.ones_like(h) * background_internal_energy
+ internal_energy[central_particle] = heated_internal_energy
-# Now we have all the information we need to set up the initial conditions!
-output = Writer(unit_system=unit_system, box_size=side_length)
+ # Now we have all the information we need to set up the initial conditions!
+ output = Writer(unit_system=unit_system, box_size=side_length)
-output.gas.coordinates = positions * side_length
-output.gas.velocities = np.zeros_like(positions) * cm / s
-output.gas.smoothing_length = h * side_length
-output.gas.internal_energy = internal_energy
-output.gas.masses = np.ones_like(h) * particle_mass
+ output.gas.coordinates = positions * side_length
+ output.gas.velocities = np.zeros_like(positions) * cm / s
+ output.gas.smoothing_length = h * side_length
+ output.gas.internal_energy = internal_energy
+ output.gas.masses = np.ones_like(h) * particle_mass
-output.write(file_name)
+ output.write(file_name)
diff --git a/examples/Cooling/FeedbackEvent_3D/plotEnergy.py b/examples/Cooling/FeedbackEvent_3D/plotEnergy.py
new file mode 100644
index 0000000000000000000000000000000000000000..6a199928b34d3d5f6ac163bfec70f9610ccafddf
--- /dev/null
+++ b/examples/Cooling/FeedbackEvent_3D/plotEnergy.py
@@ -0,0 +1,79 @@
+"""
+Plots the energy from the energy.txt file for this simulation.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from swiftsimio import load
+
+from unyt import Gyr, erg, mh, kb
+
+from makeIC import gamma, initial_density, initial_temperature, inject_temperature, mu, particle_mass
+
+try:
+ plt.style.use("mnras_durham")
+except:
+ pass
+
+
+# Snapshot for grabbing the units.
+snapshot = load("feedback_0000.hdf5")
+units = snapshot.metadata.units
+energy_units = units.mass * units.length ** 2 / (units.time ** 2)
+
+data = np.loadtxt("energy.txt").T
+
+# Assign correct units to each
+
+time = data[0] * units.time
+mass = data[1] * units.mass
+total_energy = data[2] * energy_units
+kinetic_energy = data[3] * energy_units
+thermal_energy = data[4] * energy_units
+radiative_cool = data[8] * energy_units
+
+# Now we have to figure out how much energy we actually 'injected'
+background_internal_energy = (
+ (1.0 / (mu * mh)) * (kb / (gamma - 1.0)) * initial_temperature
+)
+
+heated_internal_energy = (1.0 / (mu * mh)) * (kb / (gamma - 1.0)) * inject_temperature
+
+injected_energy = (heated_internal_energy - background_internal_energy) * particle_mass
+
+# Also want to remove the 'background' energy
+n_parts = snapshot.metadata.n_gas
+total_background_energy = background_internal_energy * n_parts * particle_mass
+
+# Now we can plot
+
+fig, ax = plt.subplots()
+
+ax.plot(
+ time.to(Gyr),
+ (kinetic_energy).to(erg),
+ label="Kinetic"
+)
+
+ax.plot(
+ time.to(Gyr),
+ (thermal_energy - total_background_energy).to(erg),
+ label="Thermal"
+)
+
+ax.plot(
+ time.to(Gyr),
+ (radiative_cool ).to(erg),
+ label="Lost to cooling"
+)
+
+ax.set_xlim(0, 0.05 * Gyr)
+
+ax.set_xlabel("Time [Gyr]")
+ax.set_ylabel("Energy [erg]")
+
+ax.legend()
+
+fig.tight_layout()
+fig.savefig("Energy.pdf")
\ No newline at end of file
diff --git a/examples/Cooling/FeedbackEvent_3D/plotEnergyAll.py b/examples/Cooling/FeedbackEvent_3D/plotEnergyAll.py
new file mode 100644
index 0000000000000000000000000000000000000000..be1ec94d93a8aad59c481a6a662a1462073baec6
--- /dev/null
+++ b/examples/Cooling/FeedbackEvent_3D/plotEnergyAll.py
@@ -0,0 +1,111 @@
+"""
+Plots the energy from the energy.txt file for this simulation.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from swiftsimio import load
+
+from unyt import Gyr, erg, mh, kb, Myr
+from scipy.interpolate import interp1d
+
+from makeIC import gamma, initial_density, initial_temperature, inject_temperature, mu, particle_mass
+
+try:
+ plt.style.use("mnras_durham")
+except:
+ pass
+
+time_to_plot = 25 * Myr
+diffusion_parameters = [0.1 * x for x in range(11)]
+plot_directory_name = "default_diffmax"
+
+kinetic_energy_at_time = []
+thermal_energy_at_time = []
+radiative_energy_at_time = []
+
+for diffusion in diffusion_parameters:
+ directory_name = f"{plot_directory_name}_{diffusion:1.1f}"
+
+# Snapshot for grabbing the units.
+ snapshot = load(f"{directory_name}/feedback_0000.hdf5")
+ units = snapshot.metadata.units
+ energy_units = units.mass * units.length ** 2 / (units.time ** 2)
+
+ data = np.loadtxt(f"{directory_name}/energy.txt").T
+
+# Assign correct units to each
+
+ time = data[0] * units.time
+ mass = data[1] * units.mass
+ total_energy = data[2] * energy_units
+ kinetic_energy = data[3] * energy_units
+ thermal_energy = data[4] * energy_units
+ radiative_cool = data[8] * energy_units
+
+# Now we have to figure out how much energy we actually 'injected'
+ background_internal_energy = (
+ (1.0 / (mu * mh)) * (kb / (gamma - 1.0)) * initial_temperature
+ )
+
+ heated_internal_energy = (1.0 / (mu * mh)) * (kb / (gamma - 1.0)) * inject_temperature
+
+ injected_energy = (heated_internal_energy - background_internal_energy) * particle_mass
+
+# Also want to remove the 'background' energy
+ n_parts = snapshot.metadata.n_gas
+ total_background_energy = background_internal_energy * n_parts * particle_mass
+
+ kinetic_energy_interpolated = interp1d(
+ time.to(Myr),
+ kinetic_energy.to(erg)
+ )
+
+ thermal_energy_interpolated = interp1d(
+ time.to(Myr),
+ (thermal_energy - total_background_energy).to(erg)
+ )
+
+ radiative_cool_interpolated = interp1d(
+ time.to(Myr),
+ radiative_cool.to(erg)
+ )
+
+ kinetic_energy_at_time.append(kinetic_energy_interpolated(time_to_plot.to(Myr)))
+ thermal_energy_at_time.append(thermal_energy_interpolated(time_to_plot.to(Myr)))
+ radiative_energy_at_time.append(radiative_cool_interpolated(time_to_plot.to(Myr)))
+
+
+# Now we can plot
+
+
+fig, ax = plt.subplots()
+
+ax.plot(
+ diffusion_parameters,
+ kinetic_energy_at_time,
+ label="Kinetic"
+)
+
+ax.plot(
+ diffusion_parameters,
+ thermal_energy_at_time,
+ label="Thermal"
+)
+
+ax.plot(
+ diffusion_parameters,
+ radiative_energy_at_time,
+ label="Lost to cooling"
+)
+
+ax.set_xlim(0, 1.0)
+
+ax.set_xlabel(r"Diffusion $\alpha_{\rm max}$")
+ax.set_ylabel(f"Energy in component at $t=${time_to_plot} [erg]")
+
+ax.legend()
+
+fig.tight_layout()
+fig.savefig("EnergyFuncDiff.pdf")
diff --git a/examples/Cooling/FeedbackEvent_3D/plotSolution.py b/examples/Cooling/FeedbackEvent_3D/plotSolution.py
index b9cce959883dbfd62c57b0f25c0b515b24308ab9..fe6f93996dafee2b6e81a70d4786374d86355f6f 100644
--- a/examples/Cooling/FeedbackEvent_3D/plotSolution.py
+++ b/examples/Cooling/FeedbackEvent_3D/plotSolution.py
@@ -29,6 +29,7 @@ from scipy import stats
from unyt import cm, s, km, kpc, Pa, msun, K, keV, mh
kPa = 1000 * Pa
+plot_radius = 7 * kpc
from swiftsimio import load
@@ -155,7 +156,7 @@ log = dict(
v_r=False, v_phi=False, u=False, S=False, P=False, rho=False, visc=False, diff=False
)
ylim = dict(
- v_r=[-8, 5], u=[3500, 5500], rho=[0.02, 0.15], visc=[0, 2.0], diff=[0, 0.25],
+ v_r=[-4, 25], u=[4750, 6000], rho=[0.09, 0.16], visc=[0, 2.0], diff=[0, 0.25],
P=[3e-18, 10e-18], S=[-0.5e60, 4e60]
)
@@ -198,7 +199,7 @@ for key, label in plot.items():
axis.set_xlabel("Radius ($r$) [kpc]", labelpad=0)
axis.set_ylabel(label, labelpad=0)
- axis.set_xlim(0.0, 0.7 * boxSize.to(kpc).value)
+ axis.set_xlim(0.0, plot_radius.to(kpc))
try:
axis.set_ylim(*ylim[key])
diff --git a/examples/Cooling/FeedbackEvent_3D/run.sh b/examples/Cooling/FeedbackEvent_3D/run.sh
index b91c515266eb5fcc1ed86c21651050613188fe9a..9111b8e8cf2cb9b0b333a1b97c04dac3975fac5c 100755
--- a/examples/Cooling/FeedbackEvent_3D/run.sh
+++ b/examples/Cooling/FeedbackEvent_3D/run.sh
@@ -5,3 +5,4 @@
# Plot the solution
python plotSolution.py 5
+python plotEnergy.py
diff --git a/examples/Cooling/FeedbackEvent_3D/runs.sh b/examples/Cooling/FeedbackEvent_3D/runs.sh
new file mode 100644
index 0000000000000000000000000000000000000000..5e3853f1e196a2bedfcb32579209fd2f614cfadc
--- /dev/null
+++ b/examples/Cooling/FeedbackEvent_3D/runs.sh
@@ -0,0 +1,31 @@
+#!/bin/bash -l
+
+#SBATCH -J SWIFTDiffusionCalibration
+#SBATCH -N 1
+#SBATCH -o swift_diffusion.out
+#SBATCH -e swift_diffusion.err
+#SBATCH -p cosma
+#SBATCH -A durham
+#SBATCH --exclusive
+
+#SBATCH -t 1:00:00
+
+
+for diffusion_alpha_max in 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0;
+do
+ mkdir default_diffmax_$diffusion_alpha_max
+
+ cd default_diffmax_$diffusion_alpha_max
+
+ ../../../swift --hydro --cooling --limiter --threads=16 --param="SPH:diffusion_alpha_max:${diffusion_alpha_max}" ../feedback.yml 2>&1 | tee output.log
+
+ cd ..
+
+ mkdir nocool_diffmax_$diffusion_alpha_max
+
+ cd nocool_diffmax_$diffusion_alpha_max
+
+ ../../../swift --hydro --temperature --limiter --threads=16 --param="SPH:diffusion_alpha_max:${diffusion_alpha_max}" ../feedback.yml 2>&1 | tee output.log
+
+ cd ..
+done
diff --git a/examples/HydroTests/BlobTest_3D/makeSliceMovie.py b/examples/HydroTests/BlobTest_3D/makeSliceMovie.py
new file mode 100644
index 0000000000000000000000000000000000000000..726c4ac01b9ae5b046a8bb3b6c1b9206764907bb
--- /dev/null
+++ b/examples/HydroTests/BlobTest_3D/makeSliceMovie.py
@@ -0,0 +1,190 @@
+"""
+Makes a movie of the output of the blob test.
+
+Josh Borrow (joshua.borrow@durham.ac.uk) 2019
+
+LGPLv3
+"""
+
+from swiftsimio import load
+from swiftsimio.visualisation import slice
+
+from p_tqdm import p_map
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from matplotlib.colors import LogNorm
+from matplotlib.animation import FuncAnimation
+
+info_frames = 15
+start_frame = 0
+end_frame = 101
+resolution = 1024
+snapshot_name = "blob"
+cmap = "Spectral_r"
+text_args = dict(color="black")
+# plot = "pressure"
+# name = "Pressure $P$"
+plot = "density"
+name = "Fluid Density $\\rho$"
+
+
+def get_image(n):
+ """
+ Gets the image for snapshot n, and also returns the associated
+ SWIFT metadata object.
+ """
+ filename = f"{snapshot_name}_{n:04d}.hdf5"
+
+ data = load(filename)
+ boxsize = data.metadata.boxsize[0].value
+
+ output = np.zeros((resolution, resolution * 4), dtype=float)
+
+ x, y, z = data.gas.coordinates.value.T
+
+ # This is an oblong box but we can only make squares!
+ for box, box_edges in enumerate([[0.0, 1.1], [0.9, 2.1], [1.9, 3.1], [2.9, 4.0]]):
+ mask = np.logical_and(x >= box_edges[0], x <= box_edges[1])
+ masked_x = x[mask] - np.float64(box)
+ masked_y = y[mask]
+ masked_z = z[mask]
+
+ hsml = data.gas.smoothing_length.value[mask]
+
+ if plot == "density":
+ mass = data.gas.masses.value[mask]
+ image = slice(
+ x=masked_y,
+ y=masked_x,
+ z=masked_z,
+ m=mass,
+ h=hsml,
+ z_slice=0.5,
+ res=resolution,
+ )
+ else:
+ quantity = getattr(data.gas, plot).value[mask]
+ # Need to divide out the particle density for non-projected density quantities
+ image = scatter(
+ x=masked_y,
+ y=masked_x,
+ z=masked_z,
+ m=quantity,
+ h=hsml,
+ z_slice=0.5,
+ res=resolution,
+ ) / scatter(
+ x=masked_y,
+ y=masked_x,
+ z=masked_z,
+ m=np.ones_like(quantity),
+ h=hsml,
+ z_slice=0.5,
+ res=resolution,
+ )
+
+ output[:, box * resolution : (box + 1) * resolution] = image
+
+ return output, data.metadata
+
+
+def get_data_dump(metadata):
+ """
+ Gets a big data dump from the SWIFT metadata
+ """
+
+ try:
+ viscosity = metadata.viscosity_info
+ except:
+ viscosity = "No info"
+
+ try:
+ diffusion = metadata.diffusion_info
+ except:
+ diffusion = "No info"
+
+ output = (
+ "$\\bf{Blob}$ $\\bf{Test}$\n\n"
+ "$\\bf{SWIFT}$\n"
+ + metadata.code_info
+ + "\n\n"
+ + "$\\bf{Compiler}$\n"
+ + metadata.compiler_info
+ + "\n\n"
+ + "$\\bf{Hydrodynamics}$\n"
+ + metadata.hydro_info
+ + "\n\n"
+ + "$\\bf{Viscosity}$\n"
+ + viscosity
+ + "\n\n"
+ + "$\\bf{Diffusion}$\n"
+ + diffusion
+ )
+
+ return output
+
+
+def time_formatter(metadata):
+ return f"$t = {metadata.t:2.2f}$"
+
+
+# Generate the frames and unpack our variables
+images, metadata = zip(*p_map(get_image, list(range(start_frame, end_frame))))
+
+# The edges are funny because of the non-periodicity.
+central_region = images[0][
+ resolution // 10 : resolution - resolution // 10,
+ resolution // 10 : resolution - resolution // 10,
+]
+norm = LogNorm(vmin=np.min(central_region), vmax=np.max(central_region), clip="black")
+
+fig, ax = plt.subplots(figsize=(8 * 4, 8), dpi=resolution // 8)
+
+fig.subplots_adjust(0, 0, 1, 1)
+ax.axis("off")
+
+# Set up the initial state
+image = ax.imshow(np.zeros_like(images[0]), norm=norm, cmap=cmap, origin="lower")
+
+description_text = ax.text(
+ 0.5,
+ 0.5,
+ get_data_dump(metadata[0]),
+ va="center",
+ ha="center",
+ **text_args,
+ transform=ax.transAxes,
+)
+
+time_text = ax.text(
+ (1 - 0.025 * 0.25),
+ 0.975,
+ time_formatter(metadata[0]),
+ **text_args,
+ va="top",
+ ha="right",
+ transform=ax.transAxes,
+)
+
+info_text = ax.text(
+ 0.025 * 0.25, 0.975, name, **text_args, va="top", ha="left", transform=ax.transAxes
+)
+
+
+def animate(n):
+ # Display just our original frames at t < 0
+ if n >= 0:
+ image.set_array(images[n])
+ description_text.set_text("")
+ time_text.set_text(time_formatter(metadata[n]))
+
+ return (image,)
+
+
+animation = FuncAnimation(
+ fig, animate, range(start_frame - info_frames, end_frame), interval=40
+)
+
+animation.save("blob_slice.mp4")
diff --git a/examples/HydroTests/KelvinHelmholtz_2D/makeMovie.py b/examples/HydroTests/KelvinHelmholtz_2D/makeMovie.py
index a52784891ab4689dcd59dc27945e573e602785f3..e40ba44dedb6e43dd25f0ce7e0b5681e61048888 100644
--- a/examples/HydroTests/KelvinHelmholtz_2D/makeMovie.py
+++ b/examples/HydroTests/KelvinHelmholtz_2D/makeMovie.py
@@ -71,7 +71,7 @@ if __name__ == "__main__":
import matplotlib.pyplot as plt
- filename = "kelvinhelmholtz"
+ filename = "kelvinHelmholtz"
dpi = 512
@@ -93,7 +93,7 @@ if __name__ == "__main__":
fig, ax = plt.subplots(1, 1, figsize=(1, 1), frameon=False)
ax.axis("off") # Remove annoying black frame.
- data_x, data_y, density = load_and_extract("kelvinhelmholtz_0000.hdf5")
+ data_x, data_y, density = load_and_extract("kelvinHelmholtz_0000.hdf5")
x = np.linspace(0, 1, dpi)
y = np.linspace(0, 1, dpi)
diff --git a/examples/HydroTests/KelvinHelmholtz_2D/makeMovieSwiftsimIO.py b/examples/HydroTests/KelvinHelmholtz_2D/makeMovieSwiftsimIO.py
new file mode 100644
index 0000000000000000000000000000000000000000..a86445e0a369bb3697793be72a3053d20f597e45
--- /dev/null
+++ b/examples/HydroTests/KelvinHelmholtz_2D/makeMovieSwiftsimIO.py
@@ -0,0 +1,102 @@
+"""
+Makes a movie of the KH 2D data.
+
+You will need to run your movie with far higher time-resolution than usual to
+get a nice movie; around 450 snapshots over 6s is required.
+
+Edit this file near the bottom with the number of snaps you have.
+
+Written by Josh Borrow (joshua.borrow@durham.ac.uk)
+"""
+
+import os
+import h5py as h5
+import numpy as np
+import scipy.interpolate as si
+
+from swiftsimio import load
+from swiftsimio.visualisation import project_gas_pixel_grid
+
+def load_and_extract(filename):
+ """
+ Load the data and extract relevant info.
+ """
+
+ return load(filename)
+
+
+def make_plot(filename, array, nx, ny, dx, dy):
+ """
+ Load the data and plop it on the grid using nearest
+ neighbour searching for finding the 'correct' value of
+ the density.
+ """
+
+ data = load_and_extract(filename)
+
+ mesh = project_gas_pixel_grid(data, nx)
+
+ array.set_array(mesh)
+
+ return array,
+
+
+def frame(n, *args):
+ """
+ Make a single frame. Requires the global variables plot and dpi.
+ """
+
+ global plot, dpi
+
+ fn = "{}_{:04d}.hdf5".format(filename, n)
+
+ return make_plot(fn, plot, dpi, dpi, (0, 1), (0, 1))
+
+
+if __name__ == "__main__":
+ import matplotlib
+ matplotlib.use("Agg")
+
+ from tqdm import tqdm
+ from matplotlib.animation import FuncAnimation
+ from scipy.stats import gaussian_kde
+
+ import matplotlib.pyplot as plt
+
+ filename = "kelvinHelmholtz"
+ dpi = 512
+
+
+ # Look for the number of files in the directory.
+ i = 0
+ while True:
+ if os.path.isfile("{}_{:04d}.hdf5".format(filename, i)):
+ i += 1
+ else:
+ break
+
+ if i > 10000:
+ raise FileNotFoundError(
+ "Could not find the snapshots in the directory")
+
+ frames = tqdm(np.arange(0, i))
+
+ # Creation of first frame
+ fig, ax = plt.subplots(1, 1, figsize=(1, 1), frameon=False)
+ ax.axis("off") # Remove annoying black frame.
+
+ data = load_and_extract("kelvinHelmholtz_0000.hdf5")
+
+
+ mesh = project_gas_pixel_grid(data, dpi)
+
+ # Global variable for set_array
+ plot = ax.imshow(mesh, extent=[0, 1, 0, 1], animated=True, interpolation="none")
+
+ anim = FuncAnimation(fig, frame, frames, interval=40, blit=False)
+
+ # Remove all whitespace
+ fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=None, hspace=None)
+
+ # Actually make the movie
+ anim.save("khmovie.mp4", dpi=dpi, bitrate=4096)
diff --git a/examples/HydroTests/KelvinHelmholtz_2D/run.sh b/examples/HydroTests/KelvinHelmholtz_2D/run.sh
index 355bf052a7ad124bcb4d88254ad780a7ffa97aba..ee3cdc4f542e3be1f28ba0deff279c612b9c49b0 100755
--- a/examples/HydroTests/KelvinHelmholtz_2D/run.sh
+++ b/examples/HydroTests/KelvinHelmholtz_2D/run.sh
@@ -10,6 +10,6 @@ fi
# Run SWIFT
../../swift --hydro --threads=4 kelvinHelmholtz.yml 2>&1 | tee output.log
+
# Plot the solution
-python plotSolution.py 6
-python makeMovie.py
+python3 makeMovieSwiftsimIO.py
diff --git a/examples/HydroTests/SodShock_BCC_3D/plotSolution.py b/examples/HydroTests/SodShock_BCC_3D/plotSolution.py
index dc1b15c3eac862365d4422f95a14ffe1713f91a6..365b679991e9a3a5bbb9e9d5108066c04e497c2f 100644
--- a/examples/HydroTests/SodShock_BCC_3D/plotSolution.py
+++ b/examples/HydroTests/SodShock_BCC_3D/plotSolution.py
@@ -35,7 +35,7 @@ from analyticSolution import analytic
snap = int(sys.argv[1])
-sim = load(f"sodshock_{snap:04d}.hdf5")
+sim = load(f"sodShock_{snap:04d}.hdf5")
# Set up plotting stuff
try:
diff --git a/examples/HydroTests/SodShock_BCC_3D/sodShock.yml b/examples/HydroTests/SodShock_BCC_3D/sodShock.yml
index 373a70bf6b6782d7bdb4ed91417a13270f6521e6..816af9c9ad620ce8617340d749b8ab3e61e53ec6 100644
--- a/examples/HydroTests/SodShock_BCC_3D/sodShock.yml
+++ b/examples/HydroTests/SodShock_BCC_3D/sodShock.yml
@@ -28,9 +28,6 @@ Statistics:
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- viscosity_alpha: 1.0
- viscosity_alpha_max: 1.0
- viscosity_alpha_min: 1.0
# Parameters related to the initial conditions
InitialConditions:
diff --git a/src/hydro/AnarchyDU/hydro.h b/src/hydro/AnarchyDU/hydro.h
index 53c66e316136b172f7610dd6f066c6e841421ced..ac4ed6c2f8176669aa7364109f95de08ddbf722c 100644
--- a/src/hydro/AnarchyDU/hydro.h
+++ b/src/hydro/AnarchyDU/hydro.h
@@ -616,6 +616,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_gradient(
struct part *restrict p) {
p->viscosity.v_sig = 2.f * p->force.soundspeed;
+ p->force.alpha_visc_max_ngb = p->viscosity.alpha;
}
/**
@@ -774,11 +775,23 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
new_diffusion_alpha += alpha_diff_dt * dt_alpha;
/* Consistency checks to ensure min < alpha < max */
- new_diffusion_alpha =
- min(new_diffusion_alpha, hydro_props->diffusion.alpha_max);
new_diffusion_alpha =
max(new_diffusion_alpha, hydro_props->diffusion.alpha_min);
+ /* Now we limit in viscous flows; remove diffusion there. If we
+ * don't do that, then we end up diffusing energy away in supernovae.
+ * This is an EAGLE-specific fix. We limit based on the maximal
+ * viscous alpha over our neighbours in an attempt to keep diffusion
+ * low near to supernovae sites. */
+
+ /* This also enforces alpha_diff < alpha_diff_max */
+
+ const float viscous_diffusion_limit =
+ hydro_props->diffusion.alpha_max *
+ (1.f - p->force.alpha_visc_max_ngb / hydro_props->viscosity.alpha_max);
+
+ new_diffusion_alpha = min(new_diffusion_alpha, viscous_diffusion_limit);
+
p->diffusion.alpha = new_diffusion_alpha;
}
diff --git a/src/hydro/AnarchyDU/hydro_iact.h b/src/hydro/AnarchyDU/hydro_iact.h
index cba945ecae8cbf2971b3d0d810cd565cb8ecc8eb..19489f7460753b15961603e68c22c039f6de27d3 100644
--- a/src/hydro/AnarchyDU/hydro_iact.h
+++ b/src/hydro/AnarchyDU/hydro_iact.h
@@ -227,6 +227,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
const float delta_u_factor = (pi->u - pj->u) * r_inv;
pi->diffusion.laplace_u += pj->mass * delta_u_factor * wi_dx / pj->rho;
pj->diffusion.laplace_u -= pi->mass * delta_u_factor * wj_dx / pi->rho;
+
+ /* Set the maximal alpha from the previous step over the neighbours
+ * (this is used to limit the diffusion in hydro_prepare_force) */
+ const float alpha_i = pi->viscosity.alpha;
+ const float alpha_j = pj->viscosity.alpha;
+ pi->force.alpha_visc_max_ngb = max(pi->force.alpha_visc_max_ngb, alpha_j);
+ pj->force.alpha_visc_max_ngb = max(pj->force.alpha_visc_max_ngb, alpha_i);
}
/**
@@ -289,6 +296,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
const float delta_u_factor = (pi->u - pj->u) * r_inv;
pi->diffusion.laplace_u += pj->mass * delta_u_factor * wi_dx / pj->rho;
+
+ /* Set the maximal alpha from the previous step over the neighbours
+ * (this is used to limit the diffusion in hydro_prepare_force) */
+ const float alpha_j = pj->viscosity.alpha;
+ pi->force.alpha_visc_max_ngb = max(pi->force.alpha_visc_max_ngb, alpha_j);
}
/**
@@ -398,9 +410,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Diffusion term */
const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
- const float v_diff =
- alpha_diff * sqrtf(0.5f * fabsf(pressurei - pressurej) / rho_ij) +
- fabsf(fac_mu * r_inv * dvdr_Hubble);
+ const float v_diff = alpha_diff * 0.5f *
+ (sqrtf(2.f * fabsf(pressurei - pressurej) / rho_ij) +
+ fabsf(fac_mu * r_inv * dvdr_Hubble));
/* wi_dx + wj_dx / 2 is F_ij */
const float diff_du_term =
v_diff * (pi->u - pj->u) * (wi_dr / rhoi + wj_dr / rhoj);
@@ -520,9 +532,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Diffusion term */
const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
- const float v_diff =
- alpha_diff * sqrtf(0.5f * fabsf(pressurei - pressurej) / rho_ij) +
- fabsf(fac_mu * r_inv * dvdr_Hubble);
+ const float v_diff = alpha_diff * 0.5f *
+ (sqrtf(2.f * fabsf(pressurei - pressurej) / rho_ij) +
+ fabsf(fac_mu * r_inv * dvdr_Hubble));
/* wi_dx + wj_dx / 2 is F_ij */
const float diff_du_term =
v_diff * (pi->u - pj->u) * (wi_dr / rhoi + wj_dr / rhoj);
diff --git a/src/hydro/AnarchyDU/hydro_part.h b/src/hydro/AnarchyDU/hydro_part.h
index c30f778b80d91d576052082c6e849fa9d1a4f38e..4b4cc187a96f5111389cde8c50f535a545e9e7f5 100644
--- a/src/hydro/AnarchyDU/hydro_part.h
+++ b/src/hydro/AnarchyDU/hydro_part.h
@@ -185,6 +185,9 @@ struct part {
/*! Balsara switch */
float balsara;
+ /*! Maximal alpha (viscosity) over neighbours */
+ float alpha_visc_max_ngb;
+
} force;
};