diff --git a/.gitignore b/.gitignore
index 798b187d4b326fb24bd4ccccabf7dc84b9bb20fd..e0e0be67879a54d27efaef9379bc545658c6d1c4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -12,6 +12,7 @@ config.h.in
config.sub
ltmain.sh
libtool
+build
src/version_string.h
swift*.tar.gz
diff --git a/doc/RTD/source/GettingStarted/running_example.rst b/doc/RTD/source/GettingStarted/running_example.rst
index 9dfbdd8c8ec98ea59892a551691aa5f230052e2e..a5614fb2a2e2068a8188d82182042feb95251a73 100644
--- a/doc/RTD/source/GettingStarted/running_example.rst
+++ b/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
diff --git a/doc/RTD/source/HydroSchemes/anarchy_sph.rst b/doc/RTD/source/HydroSchemes/anarchy_sph.rst
new file mode 100644
index 0000000000000000000000000000000000000000..8d09280039c25608d3664e1a18d9b5c28d3108b6
--- /dev/null
+++ b/doc/RTD/source/HydroSchemes/anarchy_sph.rst
@@ -0,0 +1,33 @@
+.. 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
+
diff --git a/doc/RTD/source/HydroSchemes/index.rst b/doc/RTD/source/HydroSchemes/index.rst
index 462bb7378162ff1addab3212a6901412195a3377..7331331bda7c28fad2f471dedb1335f1201d4a21 100644
--- a/doc/RTD/source/HydroSchemes/index.rst
+++ b/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
diff --git a/examples/Cooling/CoolingRedshiftDependence/.gitignore b/examples/Cooling/CoolingRedshiftDependence/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..60baa9cb833f9e075739f327f4fe68477c84a093
--- /dev/null
+++ b/examples/Cooling/CoolingRedshiftDependence/.gitignore
@@ -0,0 +1 @@
+data/*
diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml
new file mode 100644
index 0000000000000000000000000000000000000000..13ed73554936abc824ee990322cfd49464897c5b
--- /dev/null
+++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml
@@ -0,0 +1,85 @@
+# 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
diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml
new file mode 100644
index 0000000000000000000000000000000000000000..5b336bc4bf86d2596ec52e84e24da2710d96cb8a
--- /dev/null
+++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml
@@ -0,0 +1,85 @@
+# 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
diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml
new file mode 100644
index 0000000000000000000000000000000000000000..8bce6498432b3782225c7e53397868fc1356a26e
--- /dev/null
+++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml
@@ -0,0 +1,76 @@
+# 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
diff --git a/examples/Cooling/CoolingRedshiftDependence/getGlass.sh b/examples/Cooling/CoolingRedshiftDependence/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..01b4474ac21666c843b7abedfa39a76948934911
--- /dev/null
+++ b/examples/Cooling/CoolingRedshiftDependence/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/gravity_glassCube_32.hdf5
diff --git a/examples/Cooling/CoolingRedshiftDependence/makeIC.py b/examples/Cooling/CoolingRedshiftDependence/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..61b87566882f1bad9704844b13ae04773a353acd
--- /dev/null
+++ b/examples/Cooling/CoolingRedshiftDependence/makeIC.py
@@ -0,0 +1,124 @@
+"""
+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)
diff --git a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..9cde36cb05b88e60b3bf3527f3857a3774bf5dca
--- /dev/null
+++ b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py
@@ -0,0 +1,223 @@
+"""
+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)
diff --git a/examples/Cooling/CoolingRedshiftDependence/run.sh b/examples/Cooling/CoolingRedshiftDependence/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c9bf2aea7edd33886219dd4db822d930baff7d47
--- /dev/null
+++ b/examples/Cooling/CoolingRedshiftDependence/run.sh
@@ -0,0 +1,31 @@
+#!/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
diff --git a/examples/Cosmology/ComovingSodShock_3D/plotSolution.py b/examples/Cosmology/ComovingSodShock_3D/plotSolution.py
index f05a385e8620b18189d2e7abca8aebb8ae65060e..d85f4be9a49d108d133928a81ea4482fa9099792 100644
--- a/examples/Cosmology/ComovingSodShock_3D/plotSolution.py
+++ b/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)
diff --git a/examples/Cosmology/ComovingSodShock_3D/sodShock.yml b/examples/Cosmology/ComovingSodShock_3D/sodShock.yml
index 2d7a5727cbbc2cd417527ce05d7a8ea8ea05dd71..b5b2ff8fd5571c4093b043c9903398a391e51758 100644
--- a/examples/Cosmology/ComovingSodShock_3D/sodShock.yml
+++ b/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.
diff --git a/examples/HydroTests/GreshoVortex_3D/plotSolution.py b/examples/HydroTests/GreshoVortex_3D/plotSolution.py
index 8fddf148bd169310d5e69ffbfa4a6de099068c69..545440c997d9ebc3ab11562d0a7d9fa143e23ed2 100644
--- a/examples/HydroTests/GreshoVortex_3D/plotSolution.py
+++ b/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)
diff --git a/examples/HydroTests/InteractingBlastWaves_1D/interactingBlastWaves.yml b/examples/HydroTests/InteractingBlastWaves_1D/interactingBlastWaves.yml
index c4960dfa2c07b6b08cd6559b1de49f27b518bf94..aaad57dc0907cbe496ecb819d6df2a3805f7babf 100644
--- a/examples/HydroTests/InteractingBlastWaves_1D/interactingBlastWaves.yml
+++ b/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
diff --git a/examples/HydroTests/Noh_1D/noh.yml b/examples/HydroTests/Noh_1D/noh.yml
index 58e13ddda8939c8fc5fa4360a498a87f1c5b189a..3a7e328df1da1477d34d56d25121cdd4e89c03a0 100644
--- a/examples/HydroTests/Noh_1D/noh.yml
+++ b/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
+
diff --git a/examples/HydroTests/Noh_2D/noh.yml b/examples/HydroTests/Noh_2D/noh.yml
index eaf991631854e9a9781f0fcee50d996f8af949cd..76ac8eb911faca76921370ee726ee8efa00000f3 100644
--- a/examples/HydroTests/Noh_2D/noh.yml
+++ b/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
diff --git a/examples/HydroTests/SedovBlast_2D/sedov.yml b/examples/HydroTests/SedovBlast_2D/sedov.yml
index b4252581d6eb3b2932a074e7545b2d308be51865..208e4ac72207140e055b3164ceaac04f7f521e4f 100644
--- a/examples/HydroTests/SedovBlast_2D/sedov.yml
+++ b/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
diff --git a/examples/HydroTests/SedovBlast_3D/plotSolution.py b/examples/HydroTests/SedovBlast_3D/plotSolution.py
index ad34695d36f1bf8e8985b883200f17d6e38a70c9..b0f2e08441b3fa550e61602ba852228a04362fbc 100644
--- a/examples/HydroTests/SedovBlast_3D/plotSolution.py
+++ b/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)
diff --git a/examples/HydroTests/SedovBlast_3D/sedov.yml b/examples/HydroTests/SedovBlast_3D/sedov.yml
index 19e8c72538a748304ca4da076458c9ae27dc8f46..00419ba94b262a1c94c0d3cd31acf70c949b9164 100644
--- a/examples/HydroTests/SedovBlast_3D/sedov.yml
+++ b/examples/HydroTests/SedovBlast_3D/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
diff --git a/examples/HydroTests/SodShock_3D/plotSolution.py b/examples/HydroTests/SodShock_3D/plotSolution.py
index 6da7193bcd3cdfb7c22a3fc6a14f91aea5cff5f7..69b2fe4887e986156ed01e0f4177d01ccbed6035 100644
--- a/examples/HydroTests/SodShock_3D/plotSolution.py
+++ b/examples/HydroTests/SodShock_3D/plotSolution.py
@@ -84,6 +84,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
@@ -108,6 +120,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
@@ -282,13 +303,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)
diff --git a/examples/HydroTests/UniformBox_3D/makeICbig.py b/examples/HydroTests/UniformBox_3D/makeICbig.py
index bd5cf627fb535595b3abb224cbc8de50589f12cf..8841b601656d6f03684e81b4012597615485304f 100644
--- a/examples/HydroTests/UniformBox_3D/makeICbig.py
+++ b/examples/HydroTests/UniformBox_3D/makeICbig.py
@@ -45,7 +45,7 @@ internalEnergy = P / ((gamma - 1.)*rho)
n_iterations = numPart / N
remainder = numPart % N
-print "Writing", numPart, "in", n_iterations, "iterations of size", N, "and a remainder of", remainder
+print("Writing", numPart, "in", n_iterations, "iterations of size", N, "and a remainder of", remainder)
#---------------------------------------------------
@@ -55,8 +55,8 @@ file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] = [numPart % (long(1)<<32), 0, 0, 0, 0, 0]
-grp.attrs["NumPart_Total_HighWord"] = [numPart / (long(1)<<32), 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] = [numPart % (int(1)<<32), 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [numPart / (int(1)<<32), 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = 1
@@ -120,7 +120,7 @@ for n in range(n_iterations):
coords = zeros((1,3))
offset += N
- print "Done", offset,"/", numPart, "(%.1f %%)"%(100*(float)(offset)/numPart)
+ print("Done", offset,"/", numPart, "(%.1f %%)"%(100*(float)(offset)/numPart))
# And now, the remainder
v = zeros((remainder, 3))
@@ -150,7 +150,7 @@ coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
ods_x[offset:offset+remainder,:] = coords
-print "Done", offset+remainder,"/", numPart
+print("Done", offset+remainder,"/", numPart)
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index 1e50c162d02e8c326f67978a7de73fda5c6aecd8..fddc9a62067d23e13b8c2229ad5d5e798600b4e0 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -122,59 +122,65 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
/* Nothing to do here? */
if (dt == 0.) return;
- /* Current energy */
- const float u_old = hydro_get_physical_internal_energy(p, xp, cosmo);
+ /* Current energy (in internal units) */
+ const float u_old_com = hydro_get_comoving_internal_energy(p, xp);
- /* Current du_dt in physical coordinates (internal units) */
- const float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo);
+ /* Y' = RHS of the comoving equation for du/dt that will be integrated
+ forward in time using dt_therm */
+ const float hydro_du_dt_com = hydro_get_comoving_internal_energy_dt(p);
/* Calculate cooling du_dt (in cgs units) */
const double cooling_du_dt_cgs =
cooling_rate_cgs(cosmo, hydro_props, cooling, p);
/* Convert to internal units */
- float cooling_du_dt =
+ const float cooling_du_dt_physical =
cooling_du_dt_cgs * cooling->conv_factor_energy_rate_from_cgs;
- /* Add cosmological term */
- cooling_du_dt *= cosmo->a * cosmo->a;
+ /* Add cosmological term to get Y_cooling' */
+ const float cooling_du_dt = cooling_du_dt_physical * cosmo->a * cosmo->a /
+ cosmo->a_factor_internal_energy;
- float total_du_dt = hydro_du_dt + cooling_du_dt;
+ /* Y_total' */
+ float total_du_dt = hydro_du_dt_com + cooling_du_dt;
/* We now need to check that we are not going to go below any of the limits */
- /* Limit imposed by the entropy floor */
- const float A_floor = entropy_floor(p, cosmo, floor_props);
- const float rho = hydro_get_physical_density(p, cosmo);
- const float u_floor = gas_internal_energy_from_entropy(rho, A_floor);
+ /* Limit imposed by the entropy floor (comoving)
+ * (Recall entropy is the same in physical and comoving frames) */
+ const float A_floor_com = entropy_floor(p, cosmo, floor_props);
+ const float rho_com = hydro_get_comoving_density(p);
+ const float u_floor_com =
+ gas_internal_energy_from_entropy(rho_com, A_floor_com);
/* Absolute minimum */
- const float u_minimal = hydro_props->minimal_internal_energy;
+ const float u_minimal_com =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
/* Largest of both limits */
- const float u_limit = max(u_minimal, u_floor);
+ const float u_limit_com = max(u_minimal_com, u_floor_com);
/* First, check whether we may end up below the minimal energy after
* this step 1/2 kick + another 1/2 kick that could potentially be for
* a time-step twice as big. We hence check for 1.5 delta_t. */
- if (u_old + total_du_dt * 1.5 * dt_therm < u_limit) {
- total_du_dt = (u_limit - u_old) / (1.5f * dt_therm);
+ if (u_old_com + total_du_dt * 1.5 * dt_therm < u_limit_com) {
+ total_du_dt = (u_limit_com - u_old_com) / (1.5f * dt_therm);
}
/* Second, check whether the energy used in the prediction could get negative.
* We need to check for the 1/2 dt kick followed by a full time-step drift
* that could potentially be for a time-step twice as big. We hence check
* for 2.5 delta_t but this time against 0 energy not the minimum */
- if (u_old + total_du_dt * 2.5 * dt_therm < 0.) {
- total_du_dt = -u_old / ((2.5f + 0.0001f) * dt_therm);
+ if (u_old_com + total_du_dt * 2.5 * dt_therm < 0.) {
+ total_du_dt = -u_old_com / ((2.5f + 0.0001f) * dt_therm);
}
/* Update the internal energy time derivative */
- hydro_set_physical_internal_energy_dt(p, cosmo, total_du_dt);
+ hydro_set_comoving_internal_energy_dt(p, total_du_dt);
/* Store the radiated energy (assuming dt will not change) */
xp->cooling_data.radiated_energy +=
- -hydro_get_mass(p) * (total_du_dt - hydro_du_dt) * dt_therm;
+ -hydro_get_mass(p) * cooling_du_dt_physical * dt;
}
/**
@@ -208,8 +214,9 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
cooling_rate_cgs(cosmo, hydro_props, cooling, p);
/* Convert to internal units */
- const float cooling_du_dt =
- cooling_du_dt_cgs * cooling->conv_factor_energy_rate_from_cgs;
+ const float cooling_du_dt = cooling_du_dt_cgs *
+ cooling->conv_factor_energy_rate_from_cgs /
+ cosmo->a_factor_internal_energy;
/* If we are close to (or below) the limit, we ignore the condition */
if (u < 1.01f * hydro_props->minimal_internal_energy || cooling_du_dt == 0.f)
@@ -270,7 +277,7 @@ INLINE static float cooling_get_temperature(
const double mu_ionised = hydro_props->mu_ionised;
/* Particle temperature */
- const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
+ const double u = hydro_get_drifted_physical_internal_energy(p, cosmo);
/* Temperature over mean molecular weight */
const double T_over_mu = hydro_gamma_minus_one * u * m_H / k_B;
@@ -353,6 +360,9 @@ static INLINE void cooling_print_backend(
"cm^3]",
cooling->lambda_nH2_cgs);
+ message("Lambda/n_H^2=%g [internal units]",
+ cooling->lambda_nH2_cgs * cooling->conv_factor_energy_rate_from_cgs);
+
if (cooling->cooling_tstep_mult == FLT_MAX)
message("Cooling function time-step size is unlimited");
else
diff --git a/src/cooling/const_lambda/cooling_io.h b/src/cooling/const_lambda/cooling_io.h
index 9437f0f94db41725d6715cf349843bf079137305..2e2ba799ab51a73c610701499ef61f1b398e42c5 100644
--- a/src/cooling/const_lambda/cooling_io.h
+++ b/src/cooling/const_lambda/cooling_io.h
@@ -76,7 +76,10 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
UNIT_CONV_TEMPERATURE, parts,
xparts, convert_part_T);
- return 1;
+ list[1] = io_make_output_field("RadiatedEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
+ xparts, cooling_data.radiated_energy);
+
+ return 2;
}
#endif /* SWIFT_COOLING_CONST_LAMBDA_IO_H */
diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h
index 37d182889b703df9ea6d98b63a6daa5dc85274a8..0ac52165ec4b6dc5c193cf3d22d20458d2e643d3 100644
--- a/src/hydro/AnarchyPU/hydro.h
+++ b/src/hydro/AnarchyPU/hydro.h
@@ -218,6 +218,7 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) {
/* Compute the sound speed -- see theory section for justification */
/* IDEAL GAS ONLY -- P-U does not work with generic EoS. */
+
const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho);
return square_rooted;
@@ -233,7 +234,7 @@ __attribute__((always_inline)) INLINE static float
hydro_get_physical_soundspeed(const struct part *restrict p,
const struct cosmology *cosmo) {
- return cosmo->a_factor_sound_speed * p->force.soundspeed;
+ return cosmo->a_factor_sound_speed * hydro_get_comoving_soundspeed(p);
}
/**
@@ -541,8 +542,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
/* Finish calculation of the velocity divergence */
- p->viscosity.div_v *=
- h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension;
+ p->viscosity.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2;
+ p->viscosity.div_v += cosmo->H * hydro_dimension;
}
/**
@@ -693,18 +694,31 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Here we need to update the artificial viscosity */
- /* Timescale for decay */
- const float tau =
- p->h / (2.f * p->viscosity.v_sig * hydro_props->viscosity.length);
- /* Construct time differential of div.v implicitly */
- const float div_v_dt =
+ /* We use in this function that h is the radius of support */
+ const float kernel_support_physical = p->h * cosmo->a * kernel_gamma;
+ const float v_sig_physical = p->viscosity.v_sig * cosmo->a_factor_sound_speed;
+ const float soundspeed_physical = hydro_get_physical_soundspeed(p, cosmo);
+
+ const float sound_crossing_time_inverse =
+ soundspeed_physical / kernel_support_physical;
+
+ /* Construct time differential of div.v implicitly following the ANARCHY spec
+ */
+
+ float div_v_dt =
dt_alpha == 0.f
? 0.f
: (p->viscosity.div_v - p->viscosity.div_v_previous_step) / dt_alpha;
+
/* Construct the source term for the AV; if shock detected this is _positive_
* as div_v_dt should be _negative_ before the shock hits */
- const float S = p->h * p->h * max(0.f, -1.f * div_v_dt);
- const float v_sig_square = p->viscosity.v_sig * p->viscosity.v_sig;
+ const float S = kernel_support_physical * kernel_support_physical *
+ max(0.f, -1.f * div_v_dt);
+ /* 0.25 factor comes from our definition of v_sig (sum of soundspeeds rather
+ * than mean). */
+ /* Note this is v_sig_physical squared, not comoving */
+ const float v_sig_square = 0.25 * v_sig_physical * v_sig_physical;
+
/* Calculate the current appropriate value of the AV based on the above */
const float alpha_loc =
hydro_props->viscosity.alpha_max * S / (v_sig_square + S);
@@ -713,11 +727,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Reset the value of alpha to the appropriate value */
p->viscosity.alpha = alpha_loc;
} else {
- /* Integrate the alpha forward in time to decay back to alpha = 0 */
- const float alpha_dt = (alpha_loc - p->viscosity.alpha) / tau;
-
- /* Finally, we can update the actual value of the alpha */
- p->viscosity.alpha += alpha_dt * dt_alpha;
+ /* Integrate the alpha forward in time to decay back to alpha = alpha_loc */
+ p->viscosity.alpha =
+ alpha_loc + (p->viscosity.alpha - alpha_loc) *
+ expf(-dt_alpha * sound_crossing_time_inverse *
+ hydro_props->viscosity.length);
}
if (p->viscosity.alpha < hydro_props->viscosity.alpha_min) {
@@ -729,16 +743,21 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Now for the diffusive alpha */
+ const float diffusion_timescale_physical_inverse =
+ v_sig_physical / kernel_support_physical;
+
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;
+ /* Evolution term: following Schaller+ 2015. This is made up of several
+ cosmology factors: physical h, sound speed from u / sqrt(u), and the
+ 1 / a^2 coming from the laplace operator. */
+ float alpha_diff_dt = hydro_props->diffusion.beta * kernel_support_physical *
+ p->diffusion.laplace_u * cosmo->a_factor_sound_speed /
+ (sqrt_u * cosmo->a * cosmo->a);
/* 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;
+ alpha_diff_dt -= (p->diffusion.alpha - hydro_props->diffusion.alpha_min) *
+ diffusion_timescale_physical_inverse;
float new_diffusion_alpha = p->diffusion.alpha;
new_diffusion_alpha += alpha_diff_dt * dt_alpha;
diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h
index c214db3b018e00b7f3881fb301b55d6cf49a1f43..d091ebebcf4f165db006ca58667e943ddbaf8599 100644
--- a/src/hydro/AnarchyPU/hydro_iact.h
+++ b/src/hydro/AnarchyPU/hydro_iact.h
@@ -200,14 +200,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
/* We need to construct the maximal signal velocity between our particle
* and all of it's neighbours */
- const float dv_dx = (pi->v[0] - pj->v[0]) * dx[0] +
- (pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2];
+ const float r = sqrtf(r2);
+ const float r_inv = 1.f / r;
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
+
+ /* Cosmology terms for the signal velocity */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
- const float dv_dx_factor = min(0, const_viscosity_beta * dv_dx);
+ /* Add Hubble flow */
- const float new_v_sig =
- pi->force.soundspeed + pj->force.soundspeed - dv_dx_factor;
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = min(dvdr_Hubble, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Signal velocity */
+ const float new_v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Update if we need to */
pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig);
@@ -217,14 +231,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
/* Need to get some kernel values F_ij = wi_dx */
float wi, wi_dx, wj, wj_dx;
- const float r = sqrtf(r2);
const float ui = r / hi;
const float uj = r / hj;
kernel_deval(ui, &wi, &wi_dx);
kernel_deval(uj, &wj, &wj_dx);
- const float delta_u_factor = (pi->u - pj->u) / r;
+ 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;
}
@@ -253,14 +266,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
/* We need to construct the maximal signal velocity between our particle
* and all of it's neighbours */
- const float dv_dx = (pi->v[0] - pj->v[0]) * dx[0] +
- (pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2];
+ const float r = sqrtf(r2);
+ const float r_inv = 1.f / r;
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
+
+ /* Cosmology terms for the signal velocity */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
- const float dv_dx_factor = min(0, const_viscosity_beta * dv_dx);
+ /* Add Hubble flow */
- const float new_v_sig =
- pi->force.soundspeed + pj->force.soundspeed - dv_dx_factor;
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = min(dvdr_Hubble, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Signal velocity */
+ const float new_v_sig = ci + cj - const_viscosity_beta * mu_ij;
/* Update if we need to */
pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig);
@@ -269,12 +296,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
/* Need to get some kernel values F_ij = wi_dx */
float wi, wi_dx;
- const float r = sqrtf(r2);
const float ui = r / hi;
kernel_deval(ui, &wi, &wi_dx);
- const float delta_u_factor = (pi->u - pj->u) / r;
+ const float delta_u_factor = (pi->u - pj->u) * r_inv;
pi->diffusion.laplace_u += pj->mass * delta_u_factor * wi_dx / pj->rho;
}
@@ -343,7 +369,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */
- const float v_sig = 0.5 * (pi->viscosity.v_sig + pj->viscosity.v_sig);
+ const float v_sig = 0.5f * (pi->viscosity.v_sig + pj->viscosity.v_sig);
/* Balsara term */
const float balsara_i = pi->force.balsara;
@@ -390,7 +416,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Diffusion term */
const float v_diff =
max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
- const float alpha_diff = 0.5 * (pi->diffusion.alpha + pj->diffusion.alpha);
+ const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
/* wi_dx + wj_dx / 2 is F_ij */
const float diff_du_term =
alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij;
@@ -474,7 +500,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */
- const float v_sig = 0.5 * (pi->viscosity.v_sig + pj->viscosity.v_sig);
+ const float v_sig = 0.5f * (pi->viscosity.v_sig + pj->viscosity.v_sig);
/* Balsara term */
const float balsara_i = pi->force.balsara;
@@ -514,7 +540,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Diffusion term */
const float v_diff =
max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
- const float alpha_diff = 0.5 * (pi->diffusion.alpha + pj->diffusion.alpha);
+ const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
/* wi_dx + wj_dx / 2 is F_ij */
const float diff_du_term =
alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij;
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index bcef2480810c807d61a6545f04fca6144f90e8fe..f05e3229a03f18115aa0e60ba5fee94bc39c5b9e 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -562,7 +562,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->density.rot_v[2] * p->density.rot_v[2]);
/* Compute the norm of div v including the Hubble flow term */
- const float div_physical_v = p->density.div_v + 3.f * cosmo->H;
+ const float div_physical_v = p->density.div_v + hydro_dimension * cosmo->H;
const float abs_div_physical_v = fabsf(div_physical_v);
/* Compute the pressure */
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 1d5d25f734acae72490ae42552abf65983463c07..e2fd0069524de32c49cbc3cb46553b18928d14bf 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -561,7 +561,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->density.rot_v[2] * p->density.rot_v[2]);
/* Compute the norm of div v including the Hubble flow term */
- const float div_physical_v = p->density.div_v + 3.f * cosmo->H;
+ const float div_physical_v = p->density.div_v + hydro_dimension * cosmo->H;
const float abs_div_physical_v = fabsf(div_physical_v);
/* Compute the pressure */
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index 0e53b5f2ea12462fb1a618cda44934d113f655dd..3fe3c805fd4df4b495b102f9061b1e0a507e84e9 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -530,8 +530,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Finish calculation of the velocity divergence, including hubble flow term
*/
- p->density.div_v *=
- h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension;
+ p->density.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2;
+ p->density.div_v += cosmo->H * hydro_dimension;
}
/**
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
index 9f981569e80021774cebaf11066ebdac2afc6f94..be6a789ebdcd664a99a95b83f70167da04bd7534 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
@@ -526,8 +526,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
/* Finish calculation of the velocity divergence */
- p->density.div_v *=
- h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension;
+ p->density.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2;
+ p->density.div_v += cosmo->H * hydro_dimension;
}
/**
@@ -614,7 +614,14 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Artificial viscosity updates */
- const float inverse_tau = hydro_props->viscosity.length * soundspeed * h_inv;
+ /* TODO: Actually work out why this cosmology factor is correct
+ * and update the SPH / cosmology theory documents. */
+
+ /* We divide by a^2 here to make this transform under cosmology the
+ * same as the velocity (which in SWIFT has an extra 1/a^2 factor.
+ * See the cosmology theory documents for more information. */
+ const float inverse_tau =
+ (hydro_props->viscosity.length * cosmo->a2_inv) * soundspeed * h_inv;
const float source_term = -1.f * min(p->density.div_v, 0.f);
/* Compute da/dt */
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index ff578aec139ad06cb80279d8d827a41975f9773f..eb506839c859f561756fdc77072a915ff3a383aa 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -49,12 +49,15 @@
#ifdef ANARCHY_PU_SPH
/* This nasty #ifdef is only temporary until we separate the viscosity
* and hydro components. If it is not removed by July 2019, shout at JB. */
+#undef hydro_props_default_viscosity_alpha
+#define hydro_props_default_viscosity_alpha \
+ 0.1f /* Use a very low initial AV paramater for hydrodynamics tests */
#define hydro_props_default_viscosity_alpha_min \
- 0.01f /* values taken from Schaller+ 2015 */
+ 0.0f /* values NOT the same as Schaller+ 2015 */
#define hydro_props_default_viscosity_alpha_max \
2.0f /* values taken from Schaller+ 2015 */
#define hydro_props_default_viscosity_length \
- 0.01f /* values taken from Schaller+ 2015 */
+ 0.25f /* values taken from Schaller+ 2015 */
#else
#define hydro_props_default_viscosity_alpha_min \
0.1f /* values taken from (price,2004), not used in legacy gadget mode */
diff --git a/src/runner.c b/src/runner.c
index 979bdae38437129c62487e25895175db2095df70..8ac77172e1433ec0df875b0f2ed9644197f32db0 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1352,7 +1352,7 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xparts = c->hydro.xparts;
const int count = c->hydro.count;
const struct engine *e = r->e;
- const integertime_t ti_end = e->ti_current;
+ const integertime_t ti_current = e->ti_current;
const int with_cosmology = (e->policy & engine_policy_cosmology);
const double time_base = e->time_base;
const struct cosmology *cosmo = e->cosmology;
@@ -1387,9 +1387,14 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
* This is the physical time between the start and end of the time-step
* without any scale-factor powers. */
double dt_alpha;
+
if (with_cosmology) {
const integertime_t ti_step = get_integer_timestep(p->time_bin);
- dt_alpha = cosmology_get_delta_time(cosmo, ti_end - ti_step, ti_end);
+ const integertime_t ti_begin =
+ get_integer_time_begin(ti_current - 1, p->time_bin);
+
+ dt_alpha =
+ cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
} else {
dt_alpha = get_timestep(p->time_bin, time_base);
}
@@ -1568,13 +1573,16 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
* artificial viscosity and thermal conduction terms) */
const int with_cosmology = (e->policy & engine_policy_cosmology);
const double time_base = e->time_base;
- const integertime_t ti_end = e->ti_current;
+ const integertime_t ti_current = e->ti_current;
double dt_alpha;
if (with_cosmology) {
const integertime_t ti_step = get_integer_timestep(p->time_bin);
+ const integertime_t ti_begin =
+ get_integer_time_begin(ti_current - 1, p->time_bin);
+
dt_alpha =
- cosmology_get_delta_time(cosmo, ti_end - ti_step, ti_end);
+ cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
} else {
dt_alpha = get_timestep(p->time_bin, time_base);
}
@@ -1712,13 +1720,17 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
* the evolution of alpha factors (i.e. those involved in the artificial
* viscosity and thermal conduction terms) */
const int with_cosmology = (e->policy & engine_policy_cosmology);
- const integertime_t ti_end = e->ti_current;
const double time_base = e->time_base;
+ const integertime_t ti_current = e->ti_current;
double dt_alpha;
if (with_cosmology) {
const integertime_t ti_step = get_integer_timestep(p->time_bin);
- dt_alpha = cosmology_get_delta_time(cosmo, ti_end - ti_step, ti_end);
+ const integertime_t ti_begin =
+ get_integer_time_begin(ti_current - 1, p->time_bin);
+
+ dt_alpha =
+ cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
} else {
dt_alpha = get_timestep(p->time_bin, time_base);
}