diff --git a/examples/HydroTests/Diffusion_1D/.gitignore b/examples/HydroTests/Diffusion_1D/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..dadfbc6e8ce185fc82169300e6f289ed8abc9bf2
--- /dev/null
+++ b/examples/HydroTests/Diffusion_1D/.gitignore
@@ -0,0 +1 @@
+beta_*
diff --git a/examples/HydroTests/Diffusion_1D/README b/examples/HydroTests/Diffusion_1D/README
new file mode 100644
index 0000000000000000000000000000000000000000..6254e0ce9a7aa1af5f54ba388b2b08af6bb78739
--- /dev/null
+++ b/examples/HydroTests/Diffusion_1D/README
@@ -0,0 +1,19 @@
+Diffusion 1D
+============
+
+This is a very simple, 1D test. It sets up a uniform density particle
+distribution, with energy discontinuities generated at every particle.
+Particles have their internal energy set in a binary manner, with them
+alternating between "low" and "high" states. Under a standard SPH
+scheme, we expect to see no evolution whatsoever in the box. However,
+under a scheme with thermal diffusion, we expect that there will be
+some diffusion between the particles and that the distribution of
+internal energy will equalise.
+
+Included are some scripts to create initial conditions (`makeIC.py`),
+plot a solution (`plotSolution.py`), and to run the code (`run.sh`).
+Also included is a script to run a convergence test by changing the
+`SPH:beta_diffusion` parameter (`run_set.sh`).
+
+To make the initial conditions and produce the plots, the `swiftsimio`
+library (http://gitlab.cosma.dur.ac.uk/jborrow/swiftsimio) is required.
\ No newline at end of file
diff --git a/examples/HydroTests/Diffusion_1D/diffusion.yml b/examples/HydroTests/Diffusion_1D/diffusion.yml
new file mode 100644
index 0000000000000000000000000000000000000000..099cad637bf947926baa3571a2b1361d16293547
--- /dev/null
+++ b/examples/HydroTests/Diffusion_1D/diffusion.yml
@@ -0,0 +1,38 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1 # Grams
+ UnitLength_in_cgs: 1 # Centimeters
+ UnitVelocity_in_cgs: 1 # 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: 1e-0 # The end time of the simulation (in internal units).
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: diffusion # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 2e-2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-1 # 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.
+ diffusion_alpha: 0.0 # (Optional) Override the initial value for the thermal diffusion coefficient in schemes with thermal diffusion.
+ diffusion_beta: 0.01 # (Optional) Override the decay/rise rate tuning parameter for the thermal diffusion.
+ diffusion_alpha_max: 1.0 # (Optional) Override the maximal thermal diffusion coefficient that is allowed for a given particle.
+ diffusion_alpha_min: 0.0 # (Optional) Override the minimal thermal diffusion coefficient that is allowed for a given particle.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./diffusion.hdf5 # The file to read
+ periodic: 1
diff --git a/examples/HydroTests/Diffusion_1D/diffusion_fixed_alpha.yml b/examples/HydroTests/Diffusion_1D/diffusion_fixed_alpha.yml
new file mode 100644
index 0000000000000000000000000000000000000000..36f70e2492a05e90ab29496d4fc064cc2c35f62f
--- /dev/null
+++ b/examples/HydroTests/Diffusion_1D/diffusion_fixed_alpha.yml
@@ -0,0 +1,38 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1 # Grams
+ UnitLength_in_cgs: 1 # Centimeters
+ UnitVelocity_in_cgs: 1 # 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: 1e-0 # The end time of the simulation (in internal units).
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: diffusion_fixed_alpha # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 2e-2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-1 # 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.
+ diffusion_alpha: 0.1 # (Optional) Override the initial value for the thermal diffusion coefficient in schemes with thermal diffusion.
+ diffusion_beta: 0.01 # (Optional) Override the decay/rise rate tuning parameter for the thermal diffusion.
+ diffusion_alpha_max: 0.1 # (Optional) Override the maximal thermal diffusion coefficient that is allowed for a given particle.
+ diffusion_alpha_min: 0.1 # (Optional) Override the minimal thermal diffusion coefficient that is allowed for a given particle.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./diffusion.hdf5 # The file to read
+ periodic: 1
diff --git a/examples/HydroTests/Diffusion_1D/makeIC.py b/examples/HydroTests/Diffusion_1D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..09a4b5c698c7351a8a08f53a268c2b333bfef991
--- /dev/null
+++ b/examples/HydroTests/Diffusion_1D/makeIC.py
@@ -0,0 +1,143 @@
+"""
+This file is part of SWIFT.
+Copyright (c) 2019 Josh Borrow (joshua.borrow@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/>.
+
+Creates initial conditiosn for the ContactDiscontinuty_1D test.
+Requires the swiftsimio library.
+"""
+
+import numpy as np
+import unyt
+
+
+def get_particle_positions(box_length: float, n_part: int) -> np.array:
+ """
+ Gets the particle positions evenly spaced along a _periodic_ box.
+ """
+
+ dx = box_length / float(n_part)
+ x = np.arange(n_part, dtype=float) * dx + (0.5 * dx)
+
+ return x
+
+
+def get_particle_u(low: float, high: float, n_part: int) -> np.array:
+ """
+ Gets the particle internal energies, which alternate between low
+ and high. Make sure n_part is even.
+ """
+
+ indicies = np.arange(n_part)
+ u = unyt.unyt_array(np.empty(n_part, dtype=float), low.units)
+ u[...] = low
+ u[(indicies % 2).astype(bool)] = high
+
+ return u
+
+
+def get_particle_masses(mass: float, n_part: int) -> np.array:
+ """
+ Gets the particle masses.
+ """
+
+ m = unyt.unyt_array(np.empty(n_part, dtype=float), mass.units)
+ m[...] = mass
+
+ return m
+
+
+def get_particle_hsml(box_length: float, n_part: int) -> np.array:
+ """
+ Gets the particle smoothing lengths based on their MIPS.
+ """
+
+ mips = box_length / float(n_part)
+ hsml = unyt.unyt_array(np.empty(n_part, dtype=float), box_length.units)
+ hsml[...] = mips
+
+ return hsml
+
+
+if __name__ == "__main__":
+ import argparse as ap
+ from swiftsimio import Writer
+
+ parser = ap.ArgumentParser(
+ description="Makes initial conditions for the ContactDiscontinuity_1D test."
+ )
+
+ parser.add_argument(
+ "-n",
+ "--npart",
+ help="Number of particles in the box. Make sure this is even. Default: 900",
+ type=int,
+ default=900,
+ )
+
+ parser.add_argument(
+ "-m",
+ "--mass",
+ help="Particle mass in grams. Default: 1e-4",
+ type=float,
+ default=1e-4,
+ )
+
+ parser.add_argument(
+ "-l",
+ "--low",
+ help="Low value for the internal energy (in ergs/g). Default: 1e-6",
+ type=float,
+ default=1e-6,
+ )
+
+ parser.add_argument(
+ "-t",
+ "--high",
+ help="Top/high value for the internal energy (in ergs/g). Default: 1e-4",
+ type=float,
+ default=1e-4,
+ )
+
+ parser.add_argument(
+ "-b", "--boxsize", help="Boxsize in cm. Default: 1.0", type=float, default=1.0
+ )
+
+ args = vars(parser.parse_args())
+ boxsize = args["boxsize"] * unyt.cm
+ n_part = args["npart"]
+ mass = args["mass"] * unyt.g
+ low = args["low"] * unyt.erg / unyt.g
+ high = args["high"] * unyt.erg / unyt.g
+
+ if not (n_part % 2 == 0):
+ raise AttributeError("Please ensure --npart is even.")
+
+ cgs = unyt.unit_systems.cgs_unit_system
+
+ boxsize = 1.0 * unyt.cm
+
+ writer = Writer(cgs, boxsize, dimension=1)
+
+ coordinates = np.zeros((n_part, 3), dtype=float) * unyt.cm
+ coordinates[:, 0] = get_particle_positions(boxsize, n_part)
+
+ writer.gas.coordinates = coordinates
+ writer.gas.velocities = np.zeros((n_part, 3), dtype=float) * unyt.cm / unyt.s
+ writer.gas.masses = get_particle_masses(mass, n_part)
+ writer.gas.internal_energy = get_particle_u(low, high, n_part)
+ writer.gas.smoothing_length = get_particle_hsml(boxsize, n_part)
+
+ writer.write("diffusion.hdf5")
diff --git a/examples/HydroTests/Diffusion_1D/plotSolution.py b/examples/HydroTests/Diffusion_1D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..66c8ffc6418f06589a2918ae4d8ed460b0081972
--- /dev/null
+++ b/examples/HydroTests/Diffusion_1D/plotSolution.py
@@ -0,0 +1,341 @@
+"""
+This file is part of SWIFT.
+Copyright (c) 2019 Josh Borrow (joshua.borrow@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/>.
+
+Plots the solution for the ContactDiscontinuity_1D test.
+"""
+
+import matplotlib
+import matplotlib.pyplot as plt
+import numpy as np
+
+try:
+ from scipy.integrate import solve_ivp
+ solve_ode = True
+except:
+ solve_ode = False
+
+from swiftsimio import load
+
+matplotlib.use("Agg")
+
+def solve_analytic(u_0, u_1, t_0, t_1, alpha=0.1):
+ """
+ Solves the analytic equation:
+
+ $$
+ \frac{d \Delta}{d t} = \kappa \alpha (
+ \sqrt{u_0 + \Delta} + \sqrt{u_1 + \Delta}
+ ) (
+ u_1 - u_0 - 2 \Delta
+ )
+ $$
+
+ from time t0 to t1.
+
+ + alpha is the gradient term
+ + u_0 is the "low" state
+ + u_1 is the "high" state.
+ """
+
+ if not solve_ode:
+ return [0.0], [0.0]
+
+ def gradient(t, u):
+ """
+ Returns du0/dt, du1/dt
+ """
+
+ common = alpha * (np.sqrt(u[0]) + np.sqrt(u[1])) * (u[0] - u[1])
+
+ return np.array([-1.0 * common, 1.0 * common])
+
+ ret = solve_ivp(gradient, t_span=[t_0.value, t_1.value], y0=[u_0.value, u_1.value], t_eval=np.linspace(t_0.value, t_1.value, 100))
+
+ t = ret.t
+ high = ret.y[1]
+ low = ret.y[0]
+
+ return t, (high - low) * u_0.units
+
+
+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 = (
+ "SWIFT\n"
+ + metadata.code_info
+ + "\n\n"
+ + "Compiler\n"
+ + metadata.compiler_info
+ + "\n\n"
+ + "Hydrodynamics\n"
+ + metadata.hydro_info
+ + "\n\n"
+ + "Viscosity\n"
+ + viscosity
+ + "\n\n"
+ + "Diffusion\n"
+ + diffusion
+ )
+
+ return output
+
+
+def get_data_list(start: int, stop: int, handle: str):
+ """
+ Gets a list of swiftsimio objects that contains all of the data.
+ """
+
+ data = [load("{}_{:04d}.hdf5".format(handle, x)) for x in range(start, stop + 1)]
+
+ return data
+
+
+def setup_axes(size=[8, 8], dpi=300):
+ """
+ Sets up the axes with the correct labels, etc.
+ """
+ fig, ax = plt.subplots(ncols=2, nrows=2, figsize=size, dpi=dpi)
+
+ ax = ax.flatten()
+
+ ax[0].axis("off")
+ ax[1].set_xlabel("Time [s]")
+ ax[1].set_ylabel("Relative energy difference $u/\\left<u\\right>$")
+
+ ax[2].set_xlabel("Time [s]")
+ ax[2].set_ylabel("Deviation in position relative to MIPS [$\\times 10^{6}$]")
+
+ ax[3].set_xlabel("Time [s]")
+
+ return fig, ax
+
+def mean_std_max_min(data):
+ """
+ Returns:
+ mean, stdev, max, min
+ for your data.
+ """
+ means = np.array([np.mean(x) for x in data])
+ stdevs = np.array([np.std(x) for x in data])
+ maxs = np.array([np.max(x) for x in data])
+ mins = np.array([np.min(x) for x in data])
+
+ return means, stdevs, maxs, mins
+
+
+def extract_plottables_u(data_list):
+ """
+ Extracts the plottables for the internal energies. Returns:
+ mean, stdev, max, min differences between adjacent internal energies
+ """
+
+ data = [
+ np.diff(x.gas.internal_energy.value) / np.mean(x.gas.internal_energy.value)
+ for x in data_list
+ ]
+
+ return mean_std_max_min(data)
+
+
+def extract_plottables_x(data_list):
+ """
+ Extracts the plottables for positions. Returns:
+ mean, stdev, max, min * 1e6 deviations from original position
+ """
+
+ n_part = data_list[0].metadata.n_gas
+ boxsize = data_list[0].metadata.boxsize[0].value
+ dx = boxsize / n_part
+
+ original_x = np.arange(n_part, dtype=float) * dx + (0.5 * dx)
+
+ deviations = [1e6 * abs(original_x - x.gas.coordinates.value[:, 0]) / dx for x in data_list]
+
+ return mean_std_max_min(deviations)
+
+
+def extract_plottables_rho(data_list):
+ """
+ Extracts the plottables for pressure. Returns:
+ mean, stdev, max, min * 1e6 deviations from mean density
+ """
+
+ P = [x.gas.density.value for x in data_list]
+ mean_P = [np.mean(x) for x in P]
+ deviations = [1e6 * (x - y) / x for x, y in zip(mean_P, P)]
+
+ return mean_std_max_min(deviations)
+
+
+def extract_plottables_diff(data_list):
+ """
+ Extracts the plottables for pressure. Returns:
+ mean, stdev, max, min * 1e6 deviations from mean density
+ """
+
+ P = [x.gas.diffusion.value for x in data_list]
+
+ return mean_std_max_min(P)
+
+
+def make_plot(start: int, stop: int, handle: str):
+ """
+ Makes the plot and returns the figure and axes objects.
+ """
+ fig, ax = setup_axes()
+ data_list = get_data_list(start, stop, handle)
+ data_dump = get_data_dump(data_list[0].metadata)
+ t = [x.metadata.t for x in data_list]
+ means, stdevs, maxs, mins = extract_plottables_u(data_list)
+ x_means, x_stdevs, x_maxs, x_mins = extract_plottables_x(data_list)
+
+ try:
+ alpha = np.mean([np.mean(x.gas.diffusion) for x in data_list])
+ except AttributeError:
+ # Must be using a non-diffusive scheme.
+ alpha = 0.0
+
+ ax[0].text(
+ 0.5,
+ 0.5,
+ data_dump,
+ ha="center",
+ va="center",
+ fontsize=8,
+ transform=ax[0].transAxes,
+ )
+
+ ax[1].fill_between(
+ t, means - stdevs, means + stdevs, color="C0", alpha=0.5, edgecolor="none"
+ )
+ ax[1].plot(t, means, label="Mean", c="C0")
+ ax[1].plot(t, maxs, label="Max", linestyle="dashed", c="C1")
+ ax[1].plot(t, mins, label="Min", linestyle="dashed", c="C2")
+
+ if solve_ode:
+ times_ode, diff = solve_analytic(
+ u_0=data_list[0].gas.internal_energy.min(),
+ u_1=data_list[0].gas.internal_energy.max(),
+ t_0=t[0],
+ t_1=t[-1],
+ alpha=(
+ np.sqrt(5.0/3.0 * (5.0/3.0 - 1.0)) *
+ alpha / data_list[0].gas.smoothing_length[0].value
+ )
+ )
+
+ ax[1].plot(
+ times_ode,
+ (diff) / np.mean(data_list[0].gas.internal_energy),
+ label="Analytic",
+ linestyle="dotted",
+ c="C3"
+ )
+
+ #import pdb;pdb.set_trace()
+
+ ax[2].fill_between(
+ t, x_means - x_stdevs, x_means + x_stdevs, color="C0", alpha=0.5, edgecolor="none"
+ )
+ ax[2].plot(t, x_means, label="Mean", c="C0")
+ ax[2].plot(t, x_maxs, label="Max", linestyle="dashed", c="C1")
+ ax[2].plot(t, x_mins, label="Min", linestyle="dashed", c="C2")
+
+ try:
+ # Give diffusion info a go; this may not be present
+ diff_means, diff_stdevs, diff_maxs, diff_mins = extract_plottables_diff(data_list)
+
+ ax[3].set_ylabel(r"Diffusion parameter $\alpha_{diff}$")
+ ax[3].fill_between(
+ t, diff_means - diff_stdevs, diff_means + diff_stdevs, color="C0", alpha=0.5, edgecolor="none"
+ )
+ ax[3].plot(t, diff_means, label="Mean", c="C0")
+ ax[3].plot(t, diff_maxs, label="Max", linestyle="dashed", c="C1")
+ ax[3].plot(t, diff_mins, label="Min", linestyle="dashed", c="C2")
+
+ except:
+ # Diffusion info must not be present.
+ rho_means, rho_stdevs, rho_maxs, rho_mins = extract_plottables_rho(data_list)
+
+ ax[3].set_ylabel("Deviation from mean density $(\\rho_i - \\bar{\\rho}) / \\bar{\\rho}$ [$\\times 10^{6}$]")
+ ax[3].fill_between(
+ t, rho_means - rho_stdevs, rho_means + rho_stdevs, color="C0", alpha=0.5, edgecolor="none"
+ )
+ ax[3].plot(t, rho_means, label="Mean", c="C0")
+ ax[3].plot(t, rho_maxs, label="Max", linestyle="dashed", c="C1")
+ ax[3].plot(t, rho_mins, label="Min", linestyle="dashed", c="C2")
+
+ ax[1].legend(loc=1, markerfirst=False)
+
+ ax[1].set_xlim(t[0], t[-1])
+ ax[2].set_xlim(t[0], t[-1])
+ ax[3].set_xlim(t[0], t[-1])
+
+ fig.tight_layout()
+
+ return fig, ax
+
+
+if __name__ == "__main__":
+ import argparse as ap
+
+ parser = ap.ArgumentParser(
+ description="Makes a plot of the data from the ContactDiscontinuity_1D test."
+ )
+
+ parser.add_argument(
+ "-i", "--initial", help="Initial snapshot. Default: 0", type=int, default=0
+ )
+
+ parser.add_argument(
+ "-f", "--final", help="Final snapshot. Default: 50", type=int, default=50
+ )
+
+ parser.add_argument(
+ "-s",
+ "--snapshot",
+ help="First part of the snapshot filename. Default: diffusion",
+ type=str,
+ default="diffusion",
+ )
+
+ parser.add_argument(
+ "-o",
+ "--output",
+ help="Output filename. Default: diffusion.png",
+ type=str,
+ default="diffusion.png",
+ )
+
+ args = vars(parser.parse_args())
+
+ fig, ax = make_plot(args["initial"], args["final"], args["snapshot"])
+
+ fig.savefig(args["output"])
diff --git a/examples/HydroTests/Diffusion_1D/run.sh b/examples/HydroTests/Diffusion_1D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..fb70805ca01b96df20fff438b2141b03ccf6e7ab
--- /dev/null
+++ b/examples/HydroTests/Diffusion_1D/run.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e diffusion.hdf5 ]
+then
+ echo "Generating initial conditions for the Sedov blast example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../../swift --hydro --limiter --threads=1 diffusion.yml 2>&1 | tee output.log
diff --git a/examples/HydroTests/Diffusion_1D/run_set.sh b/examples/HydroTests/Diffusion_1D/run_set.sh
new file mode 100644
index 0000000000000000000000000000000000000000..d6b1f233e65b4e7851f15cf3bc8aadceff4710d7
--- /dev/null
+++ b/examples/HydroTests/Diffusion_1D/run_set.sh
@@ -0,0 +1,26 @@
+#!/bin/bash
+
+beta=(1.0 0.1 0.01 0.001)
+swift_location="../../swift"
+flags="--hydro --limiter --threads=2"
+parameter="diffusion.yml"
+parameter_fixed="diffusion_fixed_alpha.yml"
+make_ic_script="makeIC.py"
+plot_script="plotSolution.py"
+
+for i in ${beta[@]};
+do
+ mkdir beta_$i
+ cd beta_$i
+
+ python ../$make_ic_script
+
+ ../$swift_location $flags -P "SPH:diffusion_beta:${i}" ../$parameter
+ python ../$plot_script
+ ../$swift_location $flags -P "SPH:diffusion_beta:${i}" ../$parameter_fixed
+ python ../$plot_script -s diffusion_fixed_alpha -o diffusion_fixed_alpha.png
+
+ rm *.hdf5
+
+ cd ..
+done
diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h
index 9bb53f290acb16b4f9efc44e430a78a6d1f5c5ff..77c914c2b9eda678dcc4f9c41518b17938b26e1d 100644
--- a/src/hydro/AnarchyPU/hydro.h
+++ b/src/hydro/AnarchyPU/hydro.h
@@ -693,11 +693,18 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float sqrt_u = sqrtf(p->u);
/* Calculate initial value of alpha dt before bounding */
/* alpha_diff_dt is cosmology-less */
+ /* Evolution term: following Schaller+ 2015 */
float alpha_diff_dt =
hydro_props->diffusion.beta * p->h * p->diffusion.laplace_u / sqrt_u;
+ /* Decay term: not documented in Schaller+ 2015 but was present
+ * in the original EAGLE code and in Schaye+ 2015 */
+ alpha_diff_dt -=
+ (p->diffusion.alpha - hydro_props->diffusion.alpha_min) / tau;
- float new_diffusion_alpha = p->diffusion.alpha + alpha_diff_dt * dt_alpha;
+ float new_diffusion_alpha = p->diffusion.alpha;
+ new_diffusion_alpha += alpha_diff_dt * dt_alpha;
+ /* Consistency checks to ensure min < alpha < max */
if (new_diffusion_alpha > hydro_props->diffusion.alpha_max) {
new_diffusion_alpha = hydro_props->diffusion.alpha_max;
} else if (new_diffusion_alpha < hydro_props->diffusion.alpha_min) {
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index f14c88bfb5128c1da17590f50698e5d038734b71..ff578aec139ad06cb80279d8d827a41975f9773f 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -317,6 +317,12 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
io_write_attribute_f(h_grpsph, "Max v_sig ratio (limiter)",
const_limiter_max_v_sig_ratio);
+ io_write_attribute_f(h_grpsph, "Diffusion alpha", p->diffusion.alpha);
+ io_write_attribute_f(h_grpsph, "Diffusion alpha (max)",
+ p->diffusion.alpha_max);
+ io_write_attribute_f(h_grpsph, "Diffusion alpha (min)",
+ p->diffusion.alpha_min);
+ io_write_attribute_f(h_grpsph, "Diffusion beta", p->diffusion.beta);
}
#endif