Variable AV fix | Cosmology

.gitignore

@@ -12,6 +12,7 @@ config.h.in
 config.sub
 ltmain.sh
 libtool
+build
 
 src/version_string.h
 swift*.tar.gz

doc/RTD/source/GettingStarted/running_example.rst

@@ -11,7 +11,7 @@ as ``wget`` for grabbing the glass).
 
 .. code-block:: bash
    
-   cd examples/SodShock_3D
+   cd examples/HydroTests/SodShock_3D
    ./getGlass.sh
    python makeIC.py
    ../swift --hydro --threads=4 sodShock.yml

doc/RTD/source/HydroSchemes/anarchy_sph.rst

+.. ANARCHY-SPH
+   Josh Borrow 5th April 2018
+
+ANARCHY-PU SPH
+==============
+
+This scheme is similar to the one used in the EAGLE code. This scheme
+includes:
+
++ Durier & Dalla Vecchia (2012) time-step limiter
++ Pressure-Energy SPH
++ Thermal diffusion following Price (2012)
++ A simplified version of the 'Inviscid SPH' artificial viscosity
+  (Cullen & Denhen 2010).
+
+More information will be made available in a forthcoming publication.
+
+The scheme as-implemented in SWIFT is slightly different to the one
+implemented in the original EAGLE code:
+
++ Pressure-Energy SPH is used instead of Pressure-Entropy SPH
++ Artificial viscosity coefficients have changed -- from minimal
+  value of 0.1 to 0.0, and from length of 0.1 to 0.25. This
+  is based on performance of hydrodynamics tests in SWIFT and may
+  be to do with our choice of smoothing length definition.
++ Recommended kernel changed from Wendland-C2 (with 100 Ngb) to
+  Quintic Spline (with ~82 Ngb).
+
+
+.. code-block:: bash
+   
+   ./configure --with-hydro=anarchy-pu --with-kernel=quintic-spline --disable-hand-vec
+

doc/RTD/source/HydroSchemes/index.rst

@@ -17,6 +17,7 @@ schemes available in SWIFT, as well as how to implement your own.
    minimal_sph
    planetary
    hopkins_sph
+   anarchy_sph
    gizmo
    adding_your_own
 

examples/Cooling/CoolingRedshiftDependence/.gitignore

+data/*

examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun
+  UnitLength_in_cgs:   3.08567758e24 # 1 Mpc
+  UnitVelocity_in_cgs: 1e5   	     # 1 km/s
+  UnitCurrent_in_cgs:  1   	     # Amperes
+  UnitTemp_in_cgs:     1   	     # Kelvin
+
+# Cosmological parameters
+Cosmology:
+  h:              0.6777        # Reduced Hubble constant
+  a_begin:        0.5        # Initial scale-factor of the simulation (z = 1.0)
+  a_end:          0.5061559     # Final scale factor of the simulation (~ 100 myr)
+  Omega_m:        0.307         # Matter density parameter
+  Omega_lambda:   0.693         # Dark-energy density parameter
+  Omega_b:        0.0455        # Baryon density parameter
+
+# Parameters governing the time integration
+TimeIntegration:
+  dt_min:     1e-14
+  dt_max:     5e-3
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:	           data/redshift_dependence_high_z
+  delta_time:          1.000122373748838
+  scale_factor_first:  0.5
+  compression:         4
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.5
+  delta_time:          1.1
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   100      # K
+  
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./ics_high_z.hdf5     # The file to read
+  periodic:   1
+
+# Parameters for the EAGLE chemistry
+EAGLEChemistry: 	     # Solar abundances
+  init_abundance_metal:      0.
+  init_abundance_Hydrogen:   0.752
+  init_abundance_Helium:     0.248
+  init_abundance_Carbon:     0.
+  init_abundance_Nitrogen:   0.
+  init_abundance_Oxygen:     0.
+  init_abundance_Neon:       0.
+  init_abundance_Magnesium:  0.
+  init_abundance_Silicon:    0.
+  init_abundance_Iron:       0.
+
+# Parameters for the EAGLE cooling
+EAGLECooling:
+  dir_name:                ./coolingtables/
+  H_reion_z:               11.5 
+  H_reion_eV_p_H:          2.0
+  He_reion_z_centre:       3.5
+  He_reion_z_sigma:        0.5
+  He_reion_eV_p_H:         2.0
+  
+# 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
+
+LambdaCooling:
+  lambda_nH2_cgs:              2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3]
+
+ConstCooling:
+  cooling_rate: 1.          # Cooling rate (du/dt) (internal units)
+  min_energy:   1.          # Minimal internal energy per unit mass (internal units)
+  cooling_tstep_mult: 1.    # Dimensionless pre-factor for the time-step condition
\ No newline at end of file

examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e33    # 10^10 M_sun
+  UnitLength_in_cgs:   3.08567758e21 # 1 kpc
+  UnitVelocity_in_cgs: 1   	     # 1 km/s
+  UnitCurrent_in_cgs:  1   	     # Amperes
+  UnitTemp_in_cgs:     1   	     # Kelvin
+
+# Cosmological parameters
+Cosmology:
+  h:              0.6777        # Reduced Hubble constant
+  a_begin:        0.99009       # Initial scale-factor of the simulation (z = 0.01)
+  a_end:          0.99700       # Final scale factor of the simulation (~ 100 myr)
+  Omega_m:        0.307         # Matter density parameter
+  Omega_lambda:   0.693         # Dark-energy density parameter
+  Omega_b:        0.0455        # Baryon density parameter
+
+# Parameters governing the time integration
+TimeIntegration:
+  dt_min:     1e-14
+  dt_max:     5e-3
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:	           data/redshift_dependence_low_z
+  delta_time:          1.0000695206950205
+  scale_factor_first:  0.99009
+  compression:         4
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.99009
+  delta_time:          1.1
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   100      # K
+  
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./ics_low_z.hdf5     # The file to read
+  periodic:   1
+
+# Parameters for the EAGLE chemistry
+EAGLEChemistry: 	     # Primordial
+  init_abundance_metal:      0.
+  init_abundance_Hydrogen:   0.752
+  init_abundance_Helium:     0.248
+  init_abundance_Carbon:     0.
+  init_abundance_Nitrogen:   0.
+  init_abundance_Oxygen:     0.
+  init_abundance_Neon:       0.
+  init_abundance_Magnesium:  0.
+  init_abundance_Silicon:    0.
+  init_abundance_Iron:       0.
+
+# Parameters for the EAGLE cooling
+EAGLECooling:
+  dir_name:                ./coolingtables/
+  H_reion_z:               11.5 
+  H_reion_eV_p_H:          2.0
+  He_reion_z_centre:       3.5
+  He_reion_z_sigma:        0.5
+  He_reion_eV_p_H:         2.0
+  
+# 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
+
+LambdaCooling:
+  lambda_nH2_cgs:              2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3]
+
+ConstCooling:
+  cooling_rate: 1.          # Cooling rate (du/dt) (internal units)
+  min_energy:   1.          # Minimal internal energy per unit mass (internal units)
+  cooling_tstep_mult: 1.    # Dimensionless pre-factor for the time-step condition

examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun
+  UnitLength_in_cgs:   3.08567758e24 # 1 Mpc
+  UnitVelocity_in_cgs: 1e5   	     # 1 km/s
+  UnitCurrent_in_cgs:  1   	     # Amperes
+  UnitTemp_in_cgs:     1   	     # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.
+  time_end:   1.02e-4 # ~ 100 myr
+  dt_min:     1e-14
+  dt_max:     5e-5
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:	           data/redshift_dependence_no_z
+  delta_time:          1.02e-6
+  compression:         4
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1.02e-6
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   100      # K
+  
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./ics_no_z.hdf5     # The file to read
+  periodic:   1
+
+# Parameters for the EAGLE chemistry
+EAGLEChemistry: 	     # Solar abundances
+  init_abundance_metal:      0.
+  init_abundance_Hydrogen:   0.752
+  init_abundance_Helium:     0.248
+  init_abundance_Carbon:     0.
+  init_abundance_Nitrogen:   0.
+  init_abundance_Oxygen:     0.
+  init_abundance_Neon:       0.
+  init_abundance_Magnesium:  0.
+  init_abundance_Silicon:    0.
+  init_abundance_Iron:       0.
+
+# Parameters for the EAGLE cooling
+EAGLECooling:
+  dir_name:                ./coolingtables/
+  H_reion_z:               11.5 
+  H_reion_eV_p_H:          2.0
+  He_reion_z_centre:       3.5
+  He_reion_z_sigma:        0.5
+  He_reion_eV_p_H:         2.0
+  
+# 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
+
+LambdaCooling:
+  lambda_nH2_cgs:              2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3]
+
+ConstCooling:
+  cooling_rate: 1.          # Cooling rate (du/dt) (internal units)
+  min_energy:   1.          # Minimal internal energy per unit mass (internal units)
+  cooling_tstep_mult: 1.    # Dimensionless pre-factor for the time-step condition
\ No newline at end of file

examples/Cooling/CoolingRedshiftDependence/getGlass.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/gravity_glassCube_32.hdf5

examples/Cooling/CoolingRedshiftDependence/makeIC.py

+"""
+Creates a redshift-dependent cooling box such that it always has the same
+_physical_ density at the given redshift.
+"""
+
+from swiftsimio import Writer
+from swiftsimio.units import cosmo_units
+
+from unyt import mh, cm, s, K, Mpc, kb
+import numpy as np
+
+import h5py
+
+# Physics parameters.
+boxsize = 1.0 * Mpc
+physical_density = 0.1 * mh / cm ** 3
+mu_hydrogen = 0.5  # Fully ionised
+temperature = 1e7 * K
+gamma = 5.0 / 3.0
+
+
+def get_coordinates(glass_filename: str) -> np.array:
+    """
+    Gets the coordinates from the glass file.
+    """
+
+    with h5py.File(glass_filename, "r") as handle:
+        coordinates = handle["PartType1/Coordinates"][...]
+
+    return coordinates
+
+
+def generate_ics(redshift: float, filename: str, glass_filename: str) -> None:
+    """
+    Generate initial conditions for the CoolingRedshiftDependence example.
+    """
+
+    scale_factor = 1 / (1 + redshift)
+    comoving_boxsize = boxsize / scale_factor
+
+    glass_coordinates = get_coordinates(glass_filename)
+    number_of_particles = len(glass_coordinates)
+
+    gas_particle_mass = physical_density * (boxsize ** 3) / number_of_particles
+
+    writer = Writer(cosmo_units, comoving_boxsize)
+
+    writer.gas.coordinates = glass_coordinates * comoving_boxsize
+
+    writer.gas.velocities = np.zeros_like(glass_coordinates) * cm / s
+
+    writer.gas.masses = np.ones(number_of_particles, dtype=float) * gas_particle_mass
+
+    # Leave in physical units; handled by boxsize change.
+    writer.gas.internal_energy = (
+        np.ones(number_of_particles, dtype=float)
+        * 3.0 / 2.0
+        * (temperature * kb)
+        / (mu_hydrogen * mh)
+    )
+
+    writer.gas.generate_smoothing_lengths(boxsize=comoving_boxsize, dimension=3)
+
+    writer.write(filename)
+
+    return
+
+
+if __name__ == "__main__":
+    """
+    Sets up the initial parameters.
+    """
+
+    import argparse as ap
+
+    parser = ap.ArgumentParser(
+        description="""
+            Sets up the initial conditions for the cooling test. Takes two
+            redshifts, and produces two files: ics_high_z.hdf5 and
+            ics_low_z.hdf5.
+            """
+    )
+
+    parser.add_argument(
+        "-a",
+        "--high",
+        help="The high redshift to generate initial conditions for. Default: 1.0",
+        default=1,
+        type=float,
+    )
+
+    parser.add_argument(
+        "-b",
+        "--low",
+        help="The low redshift to generate initial conditions for. Default: 0.01",
+        default=0.01,
+        type=float,
+    )
+
+    parser.add_argument(
+        "-n",
+        "--nocosmo",
+        help="Generate a non-cosmological box? Default: Truthy",
+        default=1,
+        type=bool,
+    )
+
+    parser.add_argument(
+        "-g",
+        "--glass",
+        help="Glass filename. Default: gravity_glassCube_32.hdf5",
+        type=str,
+        default="gravity_glassCube_32.hdf5",
+    )
+
+    args = parser.parse_args()
+
+    generate_ics(args.low, filename="ics_low_z.hdf5", glass_filename=args.glass)
+    generate_ics(args.high, filename="ics_high_z.hdf5", glass_filename=args.glass)
+
+    if args.nocosmo:
+        generate_ics(0.0, filename="ics_no_z.hdf5", glass_filename=args.glass)
+
+    exit(0)

examples/Cooling/CoolingRedshiftDependence/plotSolution.py

+"""
+Plots the mean temperature.
+"""
+
+import matplotlib.pyplot as plt
+
+from swiftsimio import load
+
+from unyt import Myr, K, mh, cm, erg
+
+import numpy as np
+
+
+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{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 setup_axes():
+    """
+    Sets up the axes and returns fig, ax
+    """
+
+    fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(8, 8))
+
+    ax = ax.flatten()
+
+    ax[0].semilogy()
+    ax[2].semilogy()
+
+    ax[1].set_xlabel("Simulation time elapsed [Myr]")
+    ax[2].set_xlabel("Simulation time elapsed [Myr]")
+
+    ax[0].set_ylabel("Temperature of Universe [K]")
+    ax[1].set_ylabel("Physical Density of Universe [$n_H$ cm$^{-3}$]")
+    ax[2].set_ylabel("Physical Energy [erg]")
+
+    for a in ax[:-1]:
+        a.set_xlim(0, 100)
+
+    fig.tight_layout(pad=0.5)
+
+    return fig, ax
+
+
+def get_data(handle: float, n_snaps: int):
+    """
+    Returns the elapsed simulation time, temperature, and density
+    """
+
+    t0 = 0.0
+
+    times = []
+    temps = []
+    densities = []
+    energies = []
+    radiated_energies = []
+
+    for snap in range(n_snaps):
+        try:
+            data = load(f"data/{handle}_{snap:04d}.hdf5")
+
+            if snap == 0:
+                t0 = data.metadata.t.to(Myr).value
+
+            times.append(data.metadata.t.to(Myr).value - t0)
+            temps.append(np.mean(data.gas.temperature.to(K).value))
+            densities.append(
+                np.mean(data.gas.density.to(mh / cm ** 3).value)
+                / (data.metadata.scale_factor ** 3)
+            )
+
+            try:
+                energies.append(
+                    np.mean((data.gas.internal_energy * data.gas.masses).to(erg).value)
+                    * data.metadata.scale_factor ** (2)
+                )
+                radiated_energies.append(
+                    np.mean(data.gas.radiated_energy.to(erg).value)
+                )
+            except AttributeError:
+                continue
+        except OSError:
+            continue
+
+    return times, temps, densities, energies, radiated_energies
+
+
+def get_n_steps(timesteps_filename: str) -> int:
+    """
+    Reads the number of steps from the timesteps file.
+    """
+
+    data = np.genfromtxt(timesteps_filename)
+
+    return int(data[-1][0])
+
+
+def plot_single_data(
+    handle: str, n_snaps: int, label: str, ax: plt.axes, n_steps: int = 0, run: int = 0
+):
+    """
+    Takes the a single simulation and plots it on the axes.
+    """
+
+    data = get_data(handle, n_snaps)
+
+    ax[0].plot(data[0], data[1], label=label, marker="o", ms=2, c=f"C{run}")
+
+    ax[1].plot(
+        data[0], data[2], label=f"Steps: {n_steps}", marker="o", ms=2, c=f"C{run}"
+    )
+
+    if run == 0:
+        label_energy = "Particle Energy"
+        label_radiated = "Radiated Energy"
+        label_sum = "Sum"
+    else:
+        label_energy = ""
+        label_radiated = ""
+        label_sum = ""
+
+    ax[2].plot(data[0], data[3], label=label_energy, ls="dotted", C=f"C{run}")
+
+    # ax[2].plot(data[0], data[4], label=label_radiated, ls="dashed", C=f"C{run}")
+
+    # ax[2].plot(
+    #    data[0], [x + y for x, y in zip(*data[3:5])], label=label_sum, C=f"C{run}"
+    # )
+
+    return
+
+
+def make_plot(handles, names, timestep_names, n_snaps=100):
+    """
+    Makes the whole plot and returns the fig, ax objects.
+    """
+
+    fig, ax = setup_axes()
+
+    run = 0
+
+    for handle, name, timestep_name in zip(handles, names, timestep_names):
+        try:
+            n_steps = get_n_steps(timestep_name)
+        except:
+            n_steps = "Unknown"
+
+        plot_single_data(handle, n_snaps, name, ax, n_steps=n_steps, run=run)
+        run += 1
+
+    ax[0].legend()
+    ax[1].legend()
+    ax[2].legend()
+
+    info_axis = ax[-1]
+
+    try:
+        info = get_data_dump(load(f"data/{handles[0]}_0000.hdf5").metadata)
+
+        info_axis.text(
+            0.5,
+            0.45,
+            info,
+            ha="center",
+            va="center",
+            fontsize=7,
+            transform=info_axis.transAxes,
+        )
+    except OSError:
+        pass
+
+    info_axis.axis("off")
+
+    return fig, ax
+
+
+if __name__ == "__main__":
+    """
+    Plot everything!
+    """
+
+    handles = [
+        "redshift_dependence_no_z",
+        "redshift_dependence_low_z",
+        "redshift_dependence_high_z",
+    ]
+    names = ["No Cosmology", "Low Redshift ($z=0.01$)", "High Redshift ($z=1$)"]
+    timestep_names = ["timesteps_no.txt", "timesteps_low.txt", "timesteps_high.txt"]
+
+    fig, ax = make_plot(handles, names, timestep_names)
+
+    fig.savefig("redshift_dependence.png", dpi=300)

examples/Cooling/CoolingRedshiftDependence/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e gravity_glassCube_32.hdf5 ]
+then
+    echo "Fetching initial gravity glass file for the constant cosmological box example..."
+    ./getGlass.sh
+fi
+
+# Fetch the cooling tables
+if [ ! -e ics_no_z.hdf5 ]
+then
+    echo "Generating initial conditions for the uniform cosmo box example..."
+    python3 makeIC.py
+fi
+
+swift_location="../../swift"
+
+rm data/redshift_dependence_*_z_*.hdf5
+
+mkdir data
+
+# Run SWIFT
+$swift_location --hydro --cosmology --cooling --limiter --threads=4 cooling_redshift_dependence_high_z.yml 2>&1 | tee output_high.log
+mv timesteps_4.txt timesteps_high.txt
+$swift_location --hydro --cosmology --cooling --limiter --threads=4 cooling_redshift_dependence_low_z.yml 2>&1 | tee output_low.log
+mv timesteps_4.txt timesteps_low.txt
+$swift_location --hydro --cooling --limiter --threads=4 cooling_redshift_dependence_no_z.yml 2>&1 | tee output_no.log
+mv timesteps_4.txt timesteps_no.txt
+
+python3 plotSolution.py

examples/Cosmology/ComovingSodShock_3D/plotSolution.py

@@ -88,6 +88,18 @@ S = sim["/PartType0/Entropy"][:]
 P = sim["/PartType0/Pressure"][:]
 rho = sim["/PartType0/Density"][:]
 
+try:
+    diffusion = sim["/PartType0/Diffusion"][:]
+    plot_diffusion = True
+except:
+    plot_diffusion = False
+
+try:
+    viscosity = sim["/PartType0/Viscosity"][:]
+    plot_viscosity = True
+except:
+    plot_viscosity = False
+
 x_min = -1.
 x_max = 1.
 x += x_min
@@ -112,6 +124,15 @@ P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
 S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
 u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
 
+if plot_diffusion:
+    alpha_diff_bin,_,_ = stats.binned_statistic(x, diffusion, statistic='mean', bins=x_bin_edge)
+    alpha2_diff_bin,_,_ = stats.binned_statistic(x, diffusion**2, statistic='mean', bins=x_bin_edge)
+    alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin**2)
+
+if plot_viscosity:
+    alpha_visc_bin,_,_ = stats.binned_statistic(x, viscosity, statistic='mean', bins=x_bin_edge)
+    alpha2_visc_bin,_,_ = stats.binned_statistic(x, viscosity**2, statistic='mean', bins=x_bin_edge)
+    alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin**2)
 
 # Analytic solution 
 c_L = sqrt(gas_gamma * P_L / rho_L)   # Speed of the rarefaction wave
@@ -285,13 +306,26 @@ ylim(0.8, 2.2)
 
 # Entropy profile ---------------------------------
 subplot(235)
-plot(x, S, '.', color='r', ms=0.5, alpha=0.2)
-plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=0)
 xlim(-0.5, 0.5)
-ylim(0.8, 3.8)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+
+if plot_diffusion or plot_viscosity:
+    if plot_diffusion:
+        plot(x, diffusion * 100, ".", color='r', ms=0.5, alpha=0.2)
+        errorbar(x_bin, alpha_diff_bin * 100, yerr=alpha_diff_sigma_bin * 100, fmt=".", ms=8.0, color='b', lw=1.2, label="Diffusion (100x)")
+
+    if plot_viscosity:
+        plot(x, viscosity, ".", color='g', ms=0.5, alpha=0.2)
+        errorbar(x_bin, alpha_visc_bin, yerr=alpha_visc_sigma_bin, fmt=".", ms=8.0, color='y', lw=1.2, label="Viscosity")
+
+    ylabel("${\\rm{Rate~Coefficient}}~\\alpha$", labelpad=0)
+    legend()
+else:
+    plot(x, S, '.', color='r', ms=0.5, alpha=0.2)
+    plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+    errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+    ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+    ylim(0.8, 3.8)
 
 # Information -------------------------------------
 subplot(236, frameon=False)
@@ -312,5 +346,4 @@ ylim(0, 1)
 xticks([])
 yticks([])
 
-tight_layout()
 savefig("SodShock.png", dpi=200)

examples/Cosmology/ComovingSodShock_3D/sodShock.yml

@@ -2,13 +2,13 @@
 InternalUnitSystem:
   UnitMass_in_cgs:     2.94e55   # Grams
   UnitLength_in_cgs:   3.086e18   # pc
-  UnitVelocity_in_cgs: 1.   # km per s
+  UnitVelocity_in_cgs: 1   # km per s
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
 
 # Parameters governing the time integration
 TimeIntegration:
-  dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
+  dt_min:     1e-18  # 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).
 
 # Parameters governing the snapshots
@@ -16,11 +16,13 @@ Snapshots:
   basename:            sodShock # Common part of the name of output files
   time_first:          0.       # Time of the first output (in internal units)
   delta_time:          1.06638      # Time difference between consecutive outputs (in internal units)
+  #  scale_factor_first:  1.0
   scale_factor_first:  0.001
   compression:         1
   
 # Parameters governing the conserved quantities statistics
 Statistics:
+  scale_factor_first:  1.0
   delta_time:          1.02 # Time between statistics output
 
 # Parameters for the hydrodynamics scheme
@@ -41,3 +43,38 @@ Cosmology:
   a_begin: 0.001
   a_end: 0.00106638
 
+# Impose primoridal metallicity
+EAGLEChemistry:
+  init_abundance_metal:     0.
+  init_abundance_Hydrogen:  0.752
+  init_abundance_Helium:    0.248
+  init_abundance_Carbon:    0.0
+  init_abundance_Nitrogen:  0.0
+  init_abundance_Oxygen:    0.0
+  init_abundance_Neon:      0.0
+  init_abundance_Magnesium: 0.0
+  init_abundance_Silicon:   0.0
+  init_abundance_Iron:      0.0
+
+EAGLECooling:
+  dir_name:                ./coolingtables/
+  H_reion_z:               11.5 
+  H_reion_eV_p_H:          2.0
+  He_reion_z_centre:       3.5
+  He_reion_z_sigma:        0.5
+  He_reion_eV_p_H:         2.0
+
+# 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
+
+LambdaCooling:
+  lambda_nH2_cgs:              1e-48 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3])
+  cooling_tstep_mult:          1.0   # (Optional) Dimensionless pre-factor for the time-step condition.

examples/HydroTests/GreshoVortex_3D/plotSolution.py

@@ -108,6 +108,18 @@ u = sim["/PartType0/InternalEnergy"][:]
 S = sim["/PartType0/Entropy"][:]
 P = sim["/PartType0/Pressure"][:]
 
+try:
+    diffusion = sim["/PartType0/Diffusion"][:]
+    plot_diffusion = True
+except:
+    plot_diffusion = False
+
+try:
+    viscosity = sim["/PartType0/Viscosity"][:]
+    plot_viscosity = True
+except:
+    plot_viscosity = False
+
 # Bin te data
 r_bin_edge = np.arange(0., 1., 0.02)
 r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
@@ -127,6 +139,15 @@ P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
 S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
 u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
 
+if plot_diffusion:
+    alpha_diff_bin,_,_ = stats.binned_statistic(r, diffusion, statistic='mean', bins=r_bin_edge)
+    alpha2_diff_bin,_,_ = stats.binned_statistic(r, diffusion**2, statistic='mean', bins=r_bin_edge)
+    alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin**2)
+
+if plot_viscosity:
+    alpha_visc_bin,_,_ = stats.binned_statistic(r, viscosity, statistic='mean', bins=r_bin_edge)
+    alpha2_visc_bin,_,_ = stats.binned_statistic(r, viscosity**2, statistic='mean', bins=r_bin_edge)
+    alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin**2)
 
 # Plot the interesting quantities
 figure()
@@ -188,16 +209,32 @@ ylim(7.3, 9.1)
 
 # Radial entropy profile --------------------------------
 subplot(235)
-
-plot(r, S, '.', color='r', ms=0.5)
-plot(solution_r, solution_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
-plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+
 xlim(0, R_max)
-ylim(4.9 + P0, P0 + 6.1)
+
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+
+
+if plot_diffusion or plot_viscosity:
+    if plot_diffusion:
+        plot(r, diffusion, ".", color='r', ms=0.5, alpha=0.2)
+        errorbar(r_bin, alpha_diff_bin, yerr=alpha_diff_sigma_bin, fmt=".", ms=8.0, color='b', lw=1.2, label="Diffusion")
+
+    if plot_viscosity:
+        plot(r, viscosity, ".", color='g', ms=0.5, alpha=0.2)
+        errorbar(r_bin, alpha_visc_bin, yerr=alpha_visc_sigma_bin, fmt=".", ms=8.0, color='y', lw=1.2, label="Viscosity")
+
+    ylabel("${\\rm{Rate~Coefficient}}~\\alpha$", labelpad=0)
+    legend()
+else:
+    plot(r, S, '.', color='r', ms=0.5)
+    plot(solution_r, solution_s, '--', color='k', alpha=0.8, lw=1.2)
+    errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+    plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+    plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+    ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+    ylim(4.9 + P0, P0 + 6.1)
 
 # Image --------------------------------------------------
 #subplot(234)

examples/HydroTests/InteractingBlastWaves_1D/interactingBlastWaves.yml

@@ -10,7 +10,7 @@ InternalUnitSystem:
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   3.8e-2  # 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_min:     1e-10  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-5  # The maximal time-step size of the simulation (in internal units).
 
 # Parameters governing the snapshots

examples/HydroTests/Noh_1D/noh.yml

@@ -10,7 +10,7 @@ InternalUnitSystem:
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   0.6   # 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_min:     1e-10  # 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).
 
 # Parameters governing the snapshots
@@ -33,4 +33,4 @@ InitialConditions:
   file_name:  ./noh.hdf5          # The file to read
   periodic:   1
 
-  
\ No newline at end of file
+  

examples/HydroTests/Noh_2D/noh.yml

@@ -10,7 +10,7 @@ InternalUnitSystem:
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   0.6   # 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_min:     1e-10  # 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).
 
 # Parameters governing the snapshots

examples/HydroTests/SedovBlast_2D/sedov.yml

@@ -10,7 +10,7 @@ InternalUnitSystem:
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   5e-2  # 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_min:     1e-9  # 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

examples/HydroTests/SedovBlast_3D/plotSolution.py

@@ -24,7 +24,7 @@
 # Parameters
 rho_0 = 1.          # Background Density
 P_0 = 1.e-6         # Background Pressure
-E_0 = 1.            # Energy of the explosion
+E_0 = 1.0           # Energy of the explosion
 gas_gamma = 5./3.   # Gas polytropic index
 
 
@@ -86,6 +86,18 @@ S = sim["/PartType0/Entropy"][:]
 P = sim["/PartType0/Pressure"][:]
 rho = sim["/PartType0/Density"][:]
 
+try:
+    diffusion = sim["/PartType0/Diffusion"][:]
+    plot_diffusion = True
+except:
+    plot_diffusion = False
+
+try:
+    viscosity = sim["/PartType0/Viscosity"][:]
+    plot_viscosity = True
+except:
+    plot_viscosity = False
+
 # Bin te data
 r_bin_edge = np.arange(0., 0.5, 0.01)
 r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
@@ -105,6 +117,16 @@ P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
 S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
 u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
 
+if plot_diffusion:
+    alpha_diff_bin,_,_ = stats.binned_statistic(r, diffusion, statistic='mean', bins=r_bin_edge)
+    alpha2_diff_bin,_,_ = stats.binned_statistic(r, diffusion**2, statistic='mean', bins=r_bin_edge)
+    alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin**2)
+
+if plot_viscosity:
+    alpha_visc_bin,_,_ = stats.binned_statistic(r, viscosity, statistic='mean', bins=r_bin_edge)
+    alpha2_visc_bin,_,_ = stats.binned_statistic(r, viscosity**2, statistic='mean', bins=r_bin_edge)
+    alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin**2)
+
 
 # Now, work our the solution....
 
@@ -274,14 +296,28 @@ ylim(-2, 22)
 
 # Entropy profile ---------------------------------
 subplot(235)
-plot(r, S, '.', color='r', ms=0.5, alpha=0.2)
-plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=0)
-xlim(0, 1.3 * r_shock)
-ylim(-5, 50)
 
+
+if plot_diffusion or plot_viscosity:
+    if plot_diffusion:
+        plot(r, diffusion, ".", color='r', ms=0.5, alpha=0.2)
+        errorbar(r_bin, alpha_diff_bin, yerr=alpha_diff_sigma_bin, fmt=".", ms=8.0, color='b', lw=1.2, label="Diffusion")
+
+    if plot_viscosity:
+        plot(r, viscosity, ".", color='g', ms=0.5, alpha=0.2)
+        errorbar(r_bin, alpha_visc_bin, yerr=alpha_visc_sigma_bin, fmt=".", ms=8.0, color='y', lw=1.2, label="Viscosity")
+
+    ylabel("${\\rm{Rate~Coefficient}}~\\alpha$", labelpad=0)
+    legend()
+else:
+    plot(r, S, '.', color='r', ms=0.5, alpha=0.2)
+    plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+    errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+    ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+    ylim(-5, 50)
+
+xlim(0, 1.3 * r_shock)
 # Information -------------------------------------
 subplot(236, frameon=False)