Feedback Energy Injection Changes

.gitignore

@@ -30,6 +30,7 @@ examples/fof_mpi
 examples/*/*/*.xmf
 examples/*/*/*.dat
 examples/*/*/*.png
+examples/*/*/*.pdf
 examples/*/*/*.mp4
 examples/*/*/*.txt
 examples/*/*/*.rst

examples/Cooling/FeedbackEvent_3D/.gitignore

+nocool_*/*
+default_*/*

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.

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
+

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)

examples/Cooling/FeedbackEvent_3D/plotEnergy.py

+"""
+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

examples/Cooling/FeedbackEvent_3D/plotEnergyAll.py

+"""
+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")

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])

examples/Cooling/FeedbackEvent_3D/run.sh

@@ -5,3 +5,4 @@
 
 # Plot the solution
 python plotSolution.py 5
+python plotEnergy.py

examples/Cooling/FeedbackEvent_3D/runs.sh

+#!/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

examples/HydroTests/BlobTest_3D/makeSliceMovie.py

+"""
+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")

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)

examples/HydroTests/KelvinHelmholtz_2D/makeMovieSwiftsimIO.py

+"""
+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)

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

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:

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:

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;
 }
 

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);

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;
   };