Merge branch 'anarchy_fix_decay' into 'master'

Adds decay to ANARCHY and a new hydro test.

See merge request !758

examples/HydroTests/Diffusion_1D/.gitignore

+beta_*

examples/HydroTests/Diffusion_1D/README

+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

examples/HydroTests/Diffusion_1D/diffusion.yml

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

examples/HydroTests/Diffusion_1D/diffusion_fixed_alpha.yml

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

examples/HydroTests/Diffusion_1D/makeIC.py

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

examples/HydroTests/Diffusion_1D/plotSolution.py

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

examples/HydroTests/Diffusion_1D/run.sh

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

examples/HydroTests/Diffusion_1D/run_set.sh

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

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

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