Merge remote-tracking branch 'origin/master' into memreport-swift

Conflicts:
        src/engine.c

configure.ac

@@ -16,7 +16,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 # Init the project.
-AC_INIT([SWIFT],[0.8.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
+AC_INIT([SWIFT],[0.8.1],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
 swift_config_flags="$*"
 
 #  We want to stop when given unrecognised options. No subdirs so this is safe.
@@ -60,6 +60,10 @@ if test "x$ax_enable_debug" != "xno"; then
    AC_DEFINE([SWIFT_DEVELOP_MODE],1,[Enable developer code options])
 fi
 
+# C++ in GCC 6 and above has an issue with undefined the min() and max()
+# macros. This hack works around that.
+AC_DEFINE([_GLIBCXX_INCLUDE_NEXT_C_HEADERS],1,[Hack for min() and max() using g++ 6+])
+
 # Enable POSIX and platform extension preprocessor macros.
 AC_USE_SYSTEM_EXTENSIONS
 
@@ -270,6 +274,18 @@ if test "$enable_debugging_checks" = "yes"; then
    AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging])
 fi
 
+# Check if cell graph is on.
+AC_ARG_ENABLE([cell-graph],
+   [AS_HELP_STRING([--enable-cell-graph],
+     [Activate the cell graph @<:@yes/no@:>@]
+   )],
+   [enable_cell_graph="$enableval"],
+   [enable_cell_graph="no"]
+)
+if test "$enable_cell_graph" = "yes"; then
+   AC_DEFINE([SWIFT_CELL_GRAPH],1,[Enable cell graph])
+fi
+
 # Check if using our custom icbrtf is enalbled.
 AC_ARG_ENABLE([custom-icbrtf],
    [AS_HELP_STRING([--enable-custom-icbrtf],
@@ -294,6 +310,18 @@ if test "$enable_naive_interactions" = "yes"; then
    AC_DEFINE([SWIFT_USE_NAIVE_INTERACTIONS],1,[Enable use of naive cell interaction functions])
 fi
 
+# Check whether we want to default to naive cell interactions (stars)
+AC_ARG_ENABLE([naive-interactions-stars],
+   [AS_HELP_STRING([--enable-naive-interactions-stars],
+     [Activate use of naive cell interaction functions for stars @<:@yes/no@:>@]
+   )],
+   [enable_naive_interactions_stars="$enableval"],
+   [enable_naive_interactions_stars="no"]
+)
+if test "$enable_naive_interactions_stars" = "yes"; then
+   AC_DEFINE([SWIFT_USE_NAIVE_INTERACTIONS_STARS],1,[Enable use of naive cell interaction functions for stars])
+fi
+
 # Check if gravity force checks are on for some particles.
 AC_ARG_ENABLE([gravity-force-checks],
    [AS_HELP_STRING([--enable-gravity-force-checks],
@@ -1185,6 +1213,16 @@ if test "$have_armv7apmccntr"x = "yes"x; then
 	AC_DEFINE(HAVE_ARMV7A_PMCCNTR,1,[Define if you have enabled the PMCCNTR cycle counter on ARMv7a])
 fi
 
+# Check if we have native exp10 and exp10f functions. If not failback to our
+# implementations. On Apple/CLANG we have __exp10, so also check for that
+# if the compiler is clang.
+AC_CHECK_LIB([m],[exp10], [AC_DEFINE([HAVE_EXP10],1,[The exp10 function is present.])])
+AC_CHECK_LIB([m],[exp10f], [AC_DEFINE([HAVE_EXP10F],1,[The exp10f function is present.])])
+if test "$ax_cv_c_compiler_vendor" = "clang"; then
+      AC_CHECK_LIB([m],[__exp10], [AC_DEFINE([HAVE___EXP10],1,[The __exp10 function is present.])])
+      AC_CHECK_LIB([m],[__exp10f], [AC_DEFINE([HAVE___EXP10F],1,[The __exp10f function is present.])])
+fi
+
 # Add warning flags by default, if these can be used. Option =error adds
 # -Werror to GCC, clang and Intel.  Note do this last as compiler tests may
 # become errors, if that's an issue don't use CFLAGS for these, use an AC_SUBST().
@@ -1933,6 +1971,7 @@ AC_MSG_RESULT([
    Interaction debugging       : $enable_debug_interactions
    Stars interaction debugging : $enable_debug_interactions_stars
    Naive interactions          : $enable_naive_interactions
+   Naive stars interactions    : $enable_naive_interactions_stars
    Gravity checks              : $gravity_force_checks
    Custom icbrtf               : $enable_custom_icbrtf
 

doc/RTD/source/AnalysisTools/index.rst

+.. AnalysisTools
+   Loic Hausammann 20th March 2019
+
+.. _analysistools:
+
+Analysis Tools
+==============
+
+Task dependencies
+-----------------
+
+At the beginning of each simulation the file ``dependency_graph.csv`` is generated and can be transformed into a ``dot`` and a ``png`` file with the script ``tools/plot_task_dependencies.py``.
+It requires the ``dot`` package that is available in the library graphviz.
+This script has also the possibility to generate a list of function calls for each task with the option ``--with-calls`` (this list may be incomplete).
+You can convert the ``dot`` file into a ``png`` with the following command
+``dot -Tpng dependency_graph.dot -o dependency_graph.png`` or directly read it with the python module ``xdot`` with ``python -m xdot dependency_graph.dot``.
+
+
+Cell graph
+----------
+
+An interactive graph of the cells is available with the configuration option ``--enable-cell-graph``.
+During a run, SWIFT will generate a ``cell_hierarchy_*.csv`` file per MPI rank.
+The command ``tools/make_cell_hierarchy.sh cell_hierarchy_*.csv`` merges the files together and generates the file ``cell_hierarchy.html``
+that contains the graph and can be read with your favorite web browser.
+
+With chrome, you cannot access the files directly, you will need to either access them through an existing server (e.g. public http provided by your university)
+or install ``npm`` and then run the following commands
+
+.. code-block:: bash
+   
+   npm install http-server -g
+   http-server .
+
+Now you can open the web page ``http://localhost:8080/cell_hierarchy.html``.

doc/RTD/source/SubgridModels/EAGLE/index.rst

@@ -17,9 +17,10 @@ Entropy floors
 The gas particles in the EAGLE model are prevented from cooling below a
 certain temperature. The temperature limit depends on the density of the
 particles. Two floors are used in conjonction. Both are implemented as
-polytropic "equations of states" :math:`P = P_c
-\left(\rho/\rho_c\right)^\gamma`, with the constants derived from the user
-input given in terms of temperature and Hydrogen number density.
+polytropic "equations of states":math:`P = P_c
+\left(\rho/\rho_c\right)^\gamma` (all done in physical coordinates), with
+the constants derived from the user input given in terms of temperature and
+Hydrogen number density.
 
 The first limit, labelled as ``Cool``, is typically used to prevent
 low-density high-metallicity particles to cool below the warm phase because
@@ -27,12 +28,7 @@ of over-cooling induced by the absence of metal diffusion. This limit plays
 only a small role in practice. The second limit, labelled as ``Jeans``, is
 used to prevent the fragmentation of high-density gas into clumps that
 cannot be resolved by the coupled hydro+gravity solver. The two limits are
-sketched on the following figure. An additional over-density criterion is
-applied to prevent gas not collapsed into structures from being
-affected. This criterion demands that :math:`\rho > \Delta_{\rm floor}
-\Omega_b \rho_{\rm crit}`, with :math:`\Delta_{\rm floor}` specified by the
-user and :math:`\rho_{\rm crit}` the critical density at that redshift
-[#f1]_.
+sketched on the following figure.
 
 .. figure:: EAGLE_entropy_floor.svg
     :width: 400px
@@ -51,6 +47,15 @@ user and :math:`\rho_{\rm crit}` the critical density at that redshift
     the figure for clarity reasons, typical values for EAGLE runs place
     both anchors at the same temperature.
 
+An additional over-density criterion above the mean baryonic density is
+applied to prevent gas not collapsed into structures from being
+affected. To be precise, this criterion demands that the floor is applied
+only if :math:`\rho_{\rm com} > \Delta_{\rm floor}\bar{\rho_b} =
+\Delta_{\rm floor} \Omega_b \rho_{\rm crit,0}`, with :math:`\Delta_{\rm
+floor}` specified by the user, :math:`\rho_{\rm crit,0} = 3H_0/8\pi G` the
+critical density at redshift zero [#f1]_, and :math:`\rho_{\rm com}` the
+gas co-moving density. Typical values for :math:`\Delta_{\rm floor}` are of
+order 10.
 
 The model is governed by 4 parameters for each of the two
 limits. These are given in the ``EAGLEEntropyFloor`` section of the
@@ -280,7 +285,7 @@ these elements from the abundance of `Si`. More specifically, we assume that
 their abundance by mass relative to the table's solar abundance pattern is the
 same as the relative abundance of `Si` (i.e. :math:`[Ca/Si] = 0` and
 :math:`[S/Si] = 0`). Users can optionally modify the ratios used for `S` and
-`Ca`.
+`Ca`. Note that we use the *smoothed* abundances of elements for all calculations.
 
 Above the redshift of Hydrogen re-ionization we use the extra table containing
 net cooling rates for gas exposed to the CMB and a UV + X-ray background at
@@ -288,9 +293,13 @@ redshift nine truncated above 1 Rydberg. At the redshift or re-ionization, we
 additionally inject a fixed user-defined amount of energy per unit mass to all
 the gas particles.
 
-In addition to the tables we inject extra energy from Helium re-ionization using
-a Gaussian model with a user-defined redshift for the centre, width and total
-amount of energy injected per unit mass.
+In addition to the tables we inject extra energy from Helium II re-ionization
+using a Gaussian model with a user-defined redshift for the centre, width and
+total amount of energy injected per unit mass. Additional energy is also
+injected instantaneously for Hydrogen re-ionisation to all particles (active and
+inactive) to make sure the whole Universe reaches the expected temperature
+quickly (i.e not just via the interaction with the now much stronger UV
+background).
 
 For non-cosmological run, we use the :math:`z = 0` table and the interpolation
 along the redshift dimension then becomes a trivial operation.
@@ -326,7 +335,7 @@ they are listed for every gas particle:
 +---------------------+-------------------------------------+-----------+-------------------------------------+
 
 Note that if one is running without cooling switched on at runtime, the
-temperatures can be computed by passing the ``--temparature`` runtime flag (see
+temperatures can be computed by passing the ``--temperature`` runtime flag (see
 :ref:`cmdline-options`). Note that the tables then have to be available as in
 the case with cooling switched on.
 
@@ -341,9 +350,10 @@ implicit problem. A valid section of the YAML file looks like:
    EAGLECooling:
      dir_name:     /path/to/the/Wiersma/tables/directory # Absolute or relative path
      H_reion_z:            11.5      # Redhift of Hydrogen re-ionization
+     H_reion_ev_p_H:        2.0      # Energy injected in eV per Hydrogen atom for Hydrogen re-ionization.
      He_reion_z_centre:     3.5      # Centre of the Gaussian used for Helium re-ionization
      He_reion_z_sigma:      0.5      # Width of the Gaussian used for Helium re-ionization
-     He_reion_ev_p_H:       2.0      # Energy injected in eV per Hydrogen atom for Helium re-ionization.
+     He_reion_ev_p_H:       2.0      # Energy injected in eV per Hydrogen atom for Helium II re-ionization.
 
 And the optional parameters are:
 
@@ -384,8 +394,61 @@ the snapshots for each gas and star particle:
 
 .. _EAGLE_star_formation:
 
-Star formation: Schaye+2008
-~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Star formation: Schaye+2008 modified for EAGLE
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. figure:: EAGLE_SF_Z_dep.svg
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Metal-dependance of the threshold for star formation in the
+	  EAGLE model.
+
+    The dependency of the SF threshold density on the metallicty of the gas
+    in the EAGLE model (black line). The function is described by the four
+    parameters indicated on the figure. These are the slope of the
+    dependency, its position on the metallicity-axis and normalisation
+    (black circle) as well as the maximal threshold density allowed. For
+    reference, the black arrow indicates the value typically assumed for
+    solar metallicity :math:`Z_\odot=0.014` (note, however, that this value
+    does *not* enter the model at all). The values used to produce this
+    figure are the ones assumed in the reference EAGLE model.
+
+.. figure:: EAGLE_SF_EOS.svg
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Equation-of-state assumed for the star-forming gas
+
+    The equation-of-state assumed for the star-forming gas in the EAGLE
+    model (black line). The function is described by the three parameters
+    indicated on the figure. These are the slope of the relation, the
+    position of the normalisation point on the density axis and the
+    temperature expected at this density. Note that this is a normalisation
+    and *not* a threshold. Gas at densities lower than the normalisation
+    point will also be put on this equation of state when computing its
+    star formation rate. The values used to produce this figure are the
+    ones assumed in the reference EAGLE model.
+    
+.. code:: YAML
+
+   # EAGLE star formation parameters
+   EAGLEStarFormation:
+     EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
+     EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
+     EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
+     KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+     KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
+     KS_min_over_density:               57.7      # The over-density above which star-formation is allowed.
+     KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+     KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
+     KS_temperature_margin_dex:         0.5       # Logarithm base 10 of the maximal temperature difference above the EOS allowed to form stars.
+     KS_max_density_threshold_H_p_cm3:  1e5       # Hydrogen number density above which a particle gets automatically turned into a star in Hydrogen atoms per cm^3.
+     threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+     threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+     threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
+     threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+     gas_fraction:                      0.1       # The gas fraction used internally by the model.
 
 .. _EAGLE_enrichment:
 

doc/RTD/source/SubgridModels/EAGLE/plot_EAGLE_SF_EOS.py

+import matplotlib
+
+matplotlib.use("Agg")
+from pylab import *
+from scipy import stats
+
+# Plot parameters
+params = {
+    "axes.labelsize": 10,
+    "axes.titlesize": 10,
+    "font.size": 9,
+    "font.sans-serif": [
+        "Computer Modern",
+        "Computer Modern Roman",
+        "CMU Serif",
+        "cmunrm",
+        "DejaVu Sans",
+    ],
+    "mathtext.fontset": "cm",
+    "legend.fontsize": 9,
+    "xtick.labelsize": 10,
+    "ytick.labelsize": 10,
+    "text.usetex": False,
+    "figure.figsize": (3.15, 3.15),
+    "lines.markersize": 6,
+    "figure.subplot.left": 0.15,
+    "figure.subplot.right": 0.99,
+    "figure.subplot.bottom": 0.13,
+    "figure.subplot.top": 0.99,
+    "figure.subplot.wspace": 0.15,
+    "figure.subplot.hspace": 0.12,
+    "lines.linewidth": 2.0,
+}
+
+rcParams.update(params)
+
+# Equations of state
+eos_SF_rho = np.logspace(-10, 5, 1000)
+eos_SF_T = (eos_SF_rho / 10 ** (-1)) ** (1.0 / 3.0) * 8000.0
+
+# Plot the phase space diagram
+figure()
+subplot(111, xscale="log", yscale="log")
+plot(eos_SF_rho, eos_SF_T, "k-", lw=1.0)
+
+plot([1e-10, 1e-1], [8000, 8000], "k:", lw=0.6)
+plot([1e-1, 1e-1], [20, 8000], "k:", lw=0.6)
+plot([1e-1, 1e1], [20000.0, 20000.0 * 10.0 ** (2.0 / 3.0)], "k--", lw=0.6)
+text(
+    0.5e-1,
+    200000,
+    "$n_{\\rm H}$^EOS_gamma_effective",
+    va="top",
+    rotation=43,
+    fontsize=6.5,
+    family="monospace",
+)
+text(
+    0.95e-1,
+    25,
+    "EOS_density_norm_H_p_cm3",
+    rotation=90,
+    va="bottom",
+    ha="right",
+    fontsize=7,
+    family="monospace",
+)
+text(5e-8, 8400, "EOS_temperature_norm_K", va="bottom", fontsize=7, family="monospace")
+
+scatter([1e-1], [8000], s=4, color="k")
+
+xlabel("Hydrogen number density $n_{\\rm H}$ [cm$^{-3}$]", labelpad=0)
+ylabel("Temperature $T$ [K]", labelpad=2)
+xlim(3e-8, 3e3)
+ylim(20.0, 2e5)
+
+
+savefig("EAGLE_SF_EOS.svg", dpi=200)

doc/RTD/source/SubgridModels/EAGLE/plot_EAGLE_SF_Z_dep.py

+import matplotlib
+
+matplotlib.use("Agg")
+from pylab import *
+from scipy import stats
+
+# Plot parameters
+params = {
+    "axes.labelsize": 10,
+    "axes.titlesize": 10,
+    "font.size": 9,
+    "font.sans-serif": [
+        "Computer Modern",
+        "Computer Modern Roman",
+        "CMU Serif",
+        "cmunrm",
+        "DejaVu Sans",
+    ],
+    "mathtext.fontset": "cm",
+    "legend.fontsize": 9,
+    "xtick.labelsize": 10,
+    "ytick.labelsize": 10,
+    "text.usetex": False,
+    "figure.figsize": (3.15, 3.15),
+    "lines.markersize": 6,
+    "figure.subplot.left": 0.15,
+    "figure.subplot.right": 0.99,
+    "figure.subplot.bottom": 0.13,
+    "figure.subplot.top": 0.99,
+    "figure.subplot.wspace": 0.15,
+    "figure.subplot.hspace": 0.12,
+    "lines.linewidth": 2.0,
+}
+
+rcParams.update(params)
+
+# Metal dependance parameters
+Z_0 = 0.002
+norm = 0.1
+slope = -0.64
+max_n = 10.0
+
+# Function
+Z = logspace(-8, 3, 1000)
+n = norm * (Z / Z_0) ** slope
+n = np.minimum(n, np.ones(np.size(n)) * (max_n))
+
+
+# Plot the phase space diagram
+figure()
+subplot(111, xscale="log", yscale="log")
+
+plot(Z, n, "k-", lw=1.0)
+
+plot([3e-4, 3e-2], [1.0, 1.0 * 100.0 ** (slope)], "k--", lw=0.6)
+plot([1e-10, 1e10], [max_n, max_n], "k:", lw=0.6)
+plot([Z_0, Z_0], [1e-10, norm], "k:", lw=0.6)
+plot([1e-10, Z_0], [norm, norm], "k:", lw=0.6)
+scatter([Z_0], [norm], s=4, color="k")
+
+annotate(
+    "",
+    xy=(0.014, 1e-3),
+    xytext=(0.014, 3e-4),
+    arrowprops=dict(
+        facecolor="black", shrink=0.0, width=0.1, headwidth=3.0, headlength=5.0
+    ),
+)
+text(0.016, 3.5e-4, "${Z_\\odot}$", fontsize=9)
+
+text(
+    3e-4,
+    1.45,
+    "Z^threshold_slope",
+    va="top",
+    rotation=-40,
+    fontsize=7,
+    family="monospace",
+)
+text(3e-5, 12.0, "threshold_max_density_H_p_cm3", fontsize=7, family="monospace")
+text(3e-7, 0.12, "threshold_norm_H_p_cm3", fontsize=7, family="monospace")
+text(
+    0.0018,
+    0.0004,
+    "threshold_Z0",
+    rotation=90,
+    va="bottom",
+    ha="right",
+    fontsize=7,
+    family="monospace",
+)
+
+xlabel("Metallicity (metal mass fraction) $Z$ [-]", labelpad=2)
+ylabel("SF threshold number density $n_{\\rm H, thresh}$ [cm$^{-3}$]", labelpad=-1)
+
+xlim(1e-7, 1.0)
+ylim(0.0002, 50)
+
+savefig("EAGLE_SF_Z_dep.svg", dpi=200)

doc/RTD/source/SubgridModels/EAGLE/plot_EAGLE_entropy_floor.py

 import matplotlib
+
 matplotlib.use("Agg")
 from pylab import *
 from scipy import stats
@@ -8,53 +9,97 @@ params = {
     "axes.labelsize": 10,
     "axes.titlesize": 10,
     "font.size": 9,
+    "font.sans-serif": [
+        "Computer Modern",
+        "Computer Modern Roman",
+        "CMU Serif",
+        "cmunrm",
+        "DejaVu Sans",
+    ],
+    "mathtext.fontset": "cm",
     "legend.fontsize": 9,
     "xtick.labelsize": 10,
     "ytick.labelsize": 10,
-    "text.usetex": True,
+    "text.usetex": False,
     "figure.figsize": (3.15, 3.15),
+    "lines.markersize": 6,
     "figure.subplot.left": 0.15,
     "figure.subplot.right": 0.99,
     "figure.subplot.bottom": 0.13,
     "figure.subplot.top": 0.99,
     "figure.subplot.wspace": 0.15,
     "figure.subplot.hspace": 0.12,
-    "lines.markersize": 6,
     "lines.linewidth": 2.0,
-    "text.latex.unicode": True,
 }
+
 rcParams.update(params)
-rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
 
 # Equations of state
 eos_cool_rho = np.logspace(-5, 5, 1000)
-eos_cool_T = eos_cool_rho**0. * 8000.
+eos_cool_T = eos_cool_rho ** 0.0 * 8000.0
 eos_Jeans_rho = np.logspace(-1, 5, 1000)
-eos_Jeans_T = (eos_Jeans_rho/ 10**(-1))**(1./3.) * 4000.
+eos_Jeans_T = (eos_Jeans_rho / 10 ** (-1)) ** (1.0 / 3.0) * 4000.0
 
 # Plot the phase space diagram
 figure()
 subplot(111, xscale="log", yscale="log")
-plot(eos_cool_rho, eos_cool_T, 'k-', lw=1.)
-plot(eos_Jeans_rho, eos_Jeans_T, 'k-', lw=1.)
-plot([1e-10, 1e-5], [8000, 8000], 'k:', lw=0.6)
-plot([1e-10, 1e-1], [4000, 4000], 'k:', lw=0.6)
-plot([1e-5, 1e-5], [20, 8000], 'k:', lw=0.6)
-plot([1e-1, 1e-1], [20, 4000], 'k:', lw=0.6)
-plot([3e-6, 3e-4], [28000, 28000], 'k--', lw=0.6)
-text(3e-6, 22500, "$n_{\\rm H}~\\widehat{}~{\\tt Cool\\_gamma\\_effective}$", va="top", fontsize=7)
-plot([3e-1, 3e1], [15000., 15000.*10.**(2./3.)], 'k--', lw=0.6)
-text(3e-1, 200000, "$n_{\\rm H}~\\widehat{}~{\\tt Jeans\\_gamma\\_effective}$", va="top", rotation=43, fontsize=7)
-text(0.95e-5, 25, "${\\tt Cool\\_density\\_threshold\\_H\\_p\\_cm3}$", rotation=90, va="bottom", ha="right", fontsize=7)
-text(0.95e-1, 25, "${\\tt Jeans\\_density\\_threshold\\_H\\_p\\_cm3}$", rotation=90, va="bottom", ha="right", fontsize=7)
-text(5e-8, 8800, "${\\tt Cool\\_temperature\\_norm\\_K}$", va="bottom", fontsize=7)
-text(5e-8, 4400, "${\\tt Jeans\\_temperature\\_norm\\_K}$", va="bottom", fontsize=7)
-fill_between([1e-5, 1e5], [10, 10], [8000, 8000], color='0.9')
-fill_between([1e-1, 1e5], [4000, 400000], color='0.9')
-scatter([1e-5], [8000], s=4, color='k')
-scatter([1e-1], [4000], s=4, color='k')
-xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
-ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=2)
+plot(eos_cool_rho, eos_cool_T, "k-", lw=1.0)
+plot(eos_Jeans_rho, eos_Jeans_T, "k-", lw=1.0)
+plot([1e-10, 1e-5], [8000, 8000], "k:", lw=0.6)
+plot([1e-10, 1e-1], [4000, 4000], "k:", lw=0.6)
+plot([1e-5, 1e-5], [20, 8000], "k:", lw=0.6)
+plot([1e-1, 1e-1], [20, 4000], "k:", lw=0.6)
+plot([3e-6, 3e-4], [28000, 28000], "k--", lw=0.6)
+
+text(
+    3e-6,
+    22500,
+    "$n_{\\rm H}$^Cool_gamma_effective",
+    va="top",
+    fontsize=6,
+    family="monospace",
+)
+plot([3e-1, 3e1], [15000.0, 15000.0 * 10.0 ** (2.0 / 3.0)], "k--", lw=0.6)
+text(
+    3e-1,
+    200000,
+    "$n_{\\rm H}$^Jeans_gamma_effective",
+    va="top",
+    rotation=43,
+    fontsize=6,
+    family="monospace",
+)
+text(
+    0.95e-5,
+    23,
+    "Cool_density_threshold_H_p_cm3",
+    rotation=90,
+    va="bottom",
+    ha="right",
+    fontsize=6,
+    family="monospace",
+)
+text(
+    0.95e-1,
+    23,
+    "Jeans_density_threshold_H_p_cm3",
+    rotation=90,
+    va="bottom",
+    ha="right",
+    fontsize=5.5,
+    family="monospace",
+)
+text(5e-8, 8800, "Cool_temperature_norm_K", va="bottom", fontsize=6, family="monospace")
+text(
+    5e-8, 4400, "Jeans_temperature_norm_K", va="bottom", fontsize=6, family="monospace"
+)
+fill_between([1e-5, 1e5], [10, 10], [8000, 8000], color="0.9")
+fill_between([1e-1, 1e5], [4000, 400000], color="0.9")
+scatter([1e-5], [8000], s=4, color="k")
+scatter([1e-1], [4000], s=4, color="k")
+xlabel("Hydrogen number density $n_{\\rm H}$ [cm$^{-3}$]", labelpad=0)
+ylabel("Temperature $T$ [K]", labelpad=2)
 xlim(3e-8, 3e3)
-ylim(20., 2e5)
-savefig("EAGLE_entropy_floor.png", dpi=200)
+ylim(20.0, 2e5)
+
+savefig("EAGLE_entropy_floor.svg", dpi=200)

doc/RTD/source/Task/index.rst

@@ -9,11 +9,8 @@ Task System
 This section of the documentation includes information on the task system
 available in SWIFT, as well as how to implement your own task.
 
-SWIFT can produce a graph containing all the dependencies using graphviz.
-At the beginning of each simulation a ``csv`` file is generated and can be transformed into a ``png`` with the script ``tools/plot_task_dependencies.py``.
-This script has also the possibility to generate a list of function calls for each task with the option ``--with-calls``.
-You can convert the ``dot`` file into a ``png`` with the following command
-``dot -Tpng dependency_graph.dot -o dependency_graph.png`` or directly read it with the python module ``xdot`` with ``python -m xdot dependency_graph.dot``.
+SWIFT can produce a graph containing all the dependencies.
+Everything is described in :ref:`_analysistools`.
 
 
 .. toctree::

doc/RTD/source/conf.py

@@ -25,7 +25,7 @@ author = 'SWIFT Team'
 # The short X.Y version
 version = '0.8'
 # The full version, including alpha/beta/rc tags
-release = '0.8.0'
+release = '0.8.1'
 
 
 # -- General configuration ---------------------------------------------------

doc/RTD/source/index.rst

@@ -26,3 +26,4 @@ difference is the parameter file that will need to be adapted for SWIFT.
    NewOption/index
    Task/index
    VELOCIraptorInterface/index
+   AnalysisTools/index

examples/Cooling/ConstantCosmoTempEvolution/README

+This example runs a cosmological simulation using 32^3 particles with
+the gas density set to the mean baryonic density and no fluctuations
+(i.e. sigma_8 == 0). The simulation is run without gravity to prevent
+any numerical fluctuation from growing. 
+
+The gas is cooling/heating via its interaction with the UV background
+set by the cooling model. In practice, the gas will follow the
+equilibirum temperature of the model at all z and only be affected by
+reionization. Assuming primoridal abundance, the temperature of the
+gas at different redshifts can be compared to temperatures inferred
+from observations of the Lyman-alpha forrest. The plotting script runs
+this comparison once the simulation has completed.
+
+Within the EAGLE cooling model, interesting changes are to switch
+on/off the Helium II reionisation (or change its redshift) as well as
+the position and amount of energy injected by Hydrogen
+reionisation. The parameters given the YAML file are the ones used for
+actual EAGLE galaxy formation runs.

examples/Cooling/ConstantCosmoTempEvolution/const_cosmo_temp_evol.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.0099        # Initial scale-factor of the simulation (z = 100.0)
+  a_end:          1.0           # Final scale factor of the simulation
+  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-7
+  dt_max:     5e-3
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:	       cooling_box
+  delta_time:          1.04
+  scale_factor_first:  0.00991
+  compression:         4
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.00991
+  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:  ./constantBox.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

examples/Cooling/ConstantCosmoTempEvolution/datasets/schaye_et_al_2000_thermal_history.dat

+### Data from Figure 6 in Schaye et al. (2000) (2000MNRAS.318..817S)
+# z ########## z_lb ######### z_ub ######## T_0/1e4K ##### T_0_lb/1e4K ## T_0_ub/1e4K
+1.95804        1.84615        2.08042        1.09915        0.76847        2.27999
+2.20979        2.09790        2.32867        1.40112        0.90638        1.98188
+2.28322        2.22378        2.37762        1.19902        1.01236        1.36980
+2.49650        2.33217        2.62238        1.39889        1.15520        1.74657
+2.59441        2.48601        2.70629        1.60171        1.44953        1.76021
+2.65385        2.54895        2.87762        2.29650        1.57475        2.44105
+2.83916        2.70629        2.94056        1.72395        1.43950        2.30056
+3.00000        2.88462        3.13636        2.28574        1.89278        2.71473
+3.07343        2.90210        3.21329        2.01108        1.88150        2.13172
+3.23427        3.14685        3.41259        2.26386        1.40872        3.34756
+3.36014        3.22378        3.51748        1.53842        1.39604        1.69997
+3.54545        3.43007        3.69930        1.08649        0.88238        1.39077
+3.72028        3.61189        3.81469        1.25736        1.08243        1.62749
+3.83217        3.68531        4.01399        1.32092        1.03191        1.74809
+3.90909        3.81469        4.00350        1.55086        1.14616        2.01291
+4.30420        4.14685        4.42657        1.15002        0.92882        2.01412

examples/Cooling/ConstantCosmoTempEvolution/datasets/walther_et_al_2019_thermal_history.dat

+### Data from Figure 13 in Walther et al. (2019) (2019ApJ...872...13W)
+# z ######### T_0/1e4K ##### T_0_lb/1e4K ## T_0_ub/1e4K ########
+1.7793        0.77828        0.56109        1.1403 
+1.9793        0.74208        0.65158        0.90950
+2.1793        1.0226         0.87330        1.2715 
+2.3793        1.1719         0.98190        1.4570 
+2.5793        1.2398         1.0995         1.4299 
+2.7793        1.2941         1.1493         1.4796 
+2.9793        1.3032         1.1538         1.4887 
+3.1793        1.1900         1.0724         1.3258 
+3.3770        1.4118         1.2443         1.6109 
+3.5770        1.0452         0.78281        1.3575 
+3.7770        1.2081         1.0181         1.4344 
+3.9770        0.95023        0.77828        1.1674 
+4.1770        0.89593        0.82805        0.99095
+4.5770        0.88688        0.77828        1.0136 
+4.9770        0.54299        0.45701        0.66516
+5.3747        0.61086        0.47511        0.76018

examples/Cooling/ConstantCosmoTempEvolution/getGlass.sh

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

examples/Cooling/ConstantCosmoTempEvolution/makeIC.py

+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2019 Stefan Arridge (stefan.arridge@durham.ac.uk)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+from swiftsimio import Writer
+from swiftsimio.units import cosmo_units
+
+import unyt
+import numpy as np
+import h5py as h5
+
+
+# Parameters
+boxsize = 100 * unyt.Mpc
+T_i = 1000.0  # Initial temperature of the gas (in K)
+z_i = 100.0  # Initial redshift
+gamma = 5.0 / 3.0  # Gas adiabatic index
+n_p_1D = 32
+glassFile = "gravity_glassCube_32.hdf5"
+filename = "constantBox.hdf5"
+
+# Cosmology (must be same as param file)
+hubble_param = 0.6777  # same as in param file
+Omega_bar = 0.0455  # same as in param file
+
+
+# Read the glass file
+glass = h5.File(glassFile, "r")
+
+# Total number of particles
+n_p = n_p_1D ** 3
+# Read particle positions and from the glass
+glass_pos = glass["/PartType1/Coordinates"][:, :]
+glass.close()
+
+# Calculate mean baryon density today from comological parameters
+H_0 = 100.0 * hubble_param * unyt.km / unyt.s / unyt.Mpc
+rho_crit_0 = 3.0 * H_0 ** 2 / (8.0 * np.pi * unyt.G)
+rho_bar_0 = Omega_bar * rho_crit_0 
+
+# From this, we calculate the mass of the gas particles
+gas_particle_mass = rho_bar_0 * boxsize ** 3 / (n_p_1D ** 3)
+
+# Generate object. cosmo_units corresponds to default Gadget-oid units
+# of 10^10 Msun, Mpc, and km/s
+x = Writer(cosmo_units, boxsize)
+
+# 32^3 particles.
+n_p = 32 ** 3
+
+# Make gas coordinates from 0, 100 Mpc in each direction
+x.gas.coordinates = glass_pos * boxsize
+
+# Random velocities from 0 to 1 km/s
+x.gas.velocities = np.zeros((n_p, 3)) * (unyt.km / unyt.s)
+
+# Generate uniform masses as 10^6 solar masses for each particle
+x.gas.masses = np.ones(n_p, dtype=float) * gas_particle_mass
+
+# Generate internal energy corresponding to 10^4 K
+x.gas.internal_energy = (
+    np.ones(n_p, dtype=float) * (T_i * unyt.kb * unyt.K) / (1e6 * unyt.msun)
+)
+
+# Generate initial guess for smoothing lengths based on MIPS
+x.gas.generate_smoothing_lengths(boxsize=boxsize, dimension=3)
+
+# If IDs are not present, this automatically generates
+x.write(filename)

examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py

+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import h5py as h5
+import swiftsimio
+import sys
+import glob
+import unyt
+import numpy as np
+
+## read command line arguments
+snapshot_name = sys.argv[1]
+
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 9,
+'legend.fontsize': 9,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': False,
+'figure.figsize' : (4.15,3.15),
+'figure.subplot.left'    : 0.12,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.12,
+'figure.subplot.top'     : 0.99,
+'figure.subplot.wspace'  : 0.15,
+'figure.subplot.hspace'  : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 2.,
+'text.latex.unicode': True
+}
+plt.rcParams.update(params)
+
+################# Read in observational data
+
+data_schaye = np.genfromtxt("./datasets/schaye_et_al_2000_thermal_history.dat",skip_header = 2)
+data_walther = np.genfromtxt("./datasets/walther_et_al_2019_thermal_history.dat",skip_header = 2)
+
+data_schaye = data_schaye.T
+data_walther = data_walther.T
+
+schaye_z_lower_error = data_schaye[0] - data_schaye[1]
+schaye_z_upper_error = data_schaye[2] - data_schaye[0]
+schaye_T_lower_error = np.log10(data_schaye[3]*1.0e4) - np.log10(data_schaye[4]*1.0e4)
+schaye_T_upper_error = np.log10(data_schaye[5]*1.0e4) - np.log10(data_schaye[3]*1.0e4)
+walther_T_lower_error = np.log10(data_walther[1]*1.0e4) - np.log10(data_walther[2]*1.0e4)
+walther_T_upper_error = np.log10(data_walther[3]*1.0e4) - np.log10(data_walther[1]*1.0e4)
+
+############### Read in simulation data
+
+## First, get list of all snapshots
+reg_exp = "%s*.hdf5" %snapshot_name
+snap_list = glob.glob(reg_exp)
+
+z = []
+T_mean = []
+T_std = []
+rho_mean = []
+rho_std = []
+
+## loop through list
+for snap in snap_list:
+    
+    # This loads all metadata but explicitly does _not_ read any particle data yet
+    data = swiftsimio.load(snap)
+
+    # Get the redshift
+    z = np.append(z, data.metadata.z)
+        
+    # Convert gas temperatures and densities to right units
+    data.gas.temperature.convert_to_cgs()
+
+    # Get mean and standard deviation of temperature
+    T_mean.append(np.mean(data.gas.temperature) * data.gas.temperature.units)
+    T_std.append(np.std(data.gas.temperature) * data.gas.temperature.units)
+
+    # Get mean and standard deviation of density
+    rho_mean.append(np.mean(data.gas.density) * data.gas.density.units)
+    rho_std.append(np.std(data.gas.density) * data.gas.density.units)
+    
+## Turn into numpy arrays
+T_mean = np.array(T_mean) * data.gas.temperature.units
+T_std = np.array(T_std) * data.gas.temperature.units
+rho_mean = np.array(rho_mean) * data.gas.density.units
+rho_std = np.array(rho_std) * data.gas.density.units
+
+## Put Density into units of mean baryon density today
+
+# first calculate rho_bar_0 from snapshot metadata
+### from the first snapshot, get cosmology information
+d = swiftsimio.load(snap_list[0])
+cosmology = d.metadata.cosmology
+H0 = cosmology["H0 [internal units]"] / (d.units.time)
+Omega_bar = cosmology["Omega_b"]
+
+### now calculate rho_bar_0 and divide through
+rho_bar_0 = 3.0 * H0**2 / (8. * np.pi * unyt.G) * Omega_bar 
+rho_mean /= rho_bar_0
+rho_std /= rho_bar_0
+
+### sort arrays into redshift order
+ind_sorted = np.argsort(z)
+z = z[ind_sorted]
+T_mean = T_mean[ind_sorted]
+T_std = T_std[ind_sorted]
+rho_mean = rho_mean[ind_sorted]
+rho_std = rho_std[ind_sorted]
+
+### from the first snapshot, get code information
+d = swiftsimio.load(snap_list[0])
+code_info = d.metadata.code
+git_branch = code_info["Git Branch"].decode('UTF-8')
+git_revision = code_info["Git Revision"].decode('UTF-8')
+hydro_metadata = d.metadata.hydro_scheme
+scheme = hydro_metadata["Scheme"].decode('UTF-8')
+params = d.metadata.parameters
+
+subgrid_metadata = d.metadata.subgrid_scheme
+cooling_model = subgrid_metadata["Cooling Model"].decode('UTF-8')
+chemistry_model = subgrid_metadata["Chemistry Model"].decode('UTF-8')
+
+if cooling_model == 'EAGLE':
+    z_r_H = float(params['EAGLECooling:H_reion_z'])
+    H_heat_input = float(params['EAGLECooling:H_reion_eV_p_H'])
+    z_r_He_centre = float(params['EAGLECooling:He_reion_z_centre'])
+    z_r_He_sigma = float(params['EAGLECooling:He_reion_z_sigma'])
+    He_heat_input = float(params['EAGLECooling:He_reion_eV_p_H'])
+
+metallicity = "Unknown"
+if chemistry_model == 'EAGLE':
+    metallicity = float(params['EAGLEChemistry:init_abundance_metal'])
+    
+# Make plot of temperature evolution  --------------------------------
+fig = plt.figure()
+
+# Plot sim properties
+if cooling_model == 'EAGLE':
+    plt.plot([z_r_H, z_r_H], [3.4, 4.4], 'k--', alpha=0.5, lw=0.7)
+    plt.text(z_r_H + 0.1, 3.55, "H reion.", rotation=90, alpha=0.5, fontsize=7, va="bottom")
+    plt.plot([z_r_He_centre, z_r_He_centre], [3.4, 4.4], 'k--', alpha=0.5, lw=0.7)
+    plt.text(z_r_He_centre + 0.1, 3.55, "HeII reion.", rotation=90, alpha=0.5, fontsize=7, va="bottom")
+    
+# Plot observational data
+plt.errorbar(data_schaye[0],
+             np.log10(data_schaye[3]*1.0e4),
+             xerr = [schaye_z_lower_error,schaye_z_upper_error],
+             yerr = [schaye_T_lower_error,schaye_T_upper_error],
+             fmt='s', mec='0.3', color='0.3', markersize=4, markeredgewidth=0.5, linewidth=0.5, mfc='w', label="Schaye et al. (2000)")
+plt.errorbar(data_walther[0],
+             np.log10(data_walther[1]*1.0e4),
+             yerr = [walther_T_lower_error,walther_T_upper_error],
+              fmt='.', mec='0.3', color='0.3', markersize=7, markeredgewidth=0.5, linewidth=0.5, mfc='w', label = "Walther et al. (2019)")
+
+# Plot simulation
+plt.plot(z, np.log10(T_mean))
+
+# Legend
+plt.legend(loc="upper right", frameon=True, fontsize=8, handletextpad=0.1, facecolor="w", edgecolor="w", framealpha=1.)
+plt.text(0.2, 4.8, "SWIFT %s \nCooling model: %s \n$\\rho=%.3f$ $\\Omega_{b}\\rho_{crit,0}$\n$Z=%s$"%(git_revision, cooling_model, rho_mean[-1], metallicity), va="top", ha="left", fontsize=8)
+
+plt.xlim(0, 12.2)
+plt.ylim(3.5,4.85)
+plt.xlabel("Redshift", labelpad=-0.5)
+plt.ylabel(r"$\log_{10}(T/K)$", labelpad=0)
+plt.savefig("Temperature_evolution.png", dpi=200)
+
+
+# Make plot of denisty evolution  --------------------------------
+plt.rcParams.update({'figure.subplot.left'    : 0.14})
+fig = plt.figure()
+
+plt.text(0.2, 1.011, "SWIFT %s"%git_revision, va="top", ha="left", fontsize=8)
+
+plt.fill_between(z,rho_mean - rho_std,rho_mean + rho_std,alpha = 0.5)
+plt.plot(z,rho_mean)
+
+plt.axhline(y = 1.0, linestyle = '--', color='k', alpha=0.5, lw=1.)
+
+plt.xlim(0.0,12.2)
+plt.ylabel(r"$\delta_b = \rho / \Omega_b\rho_{crit,0}$", labelpad=0.)
+plt.ylim(0.988,1.012)
+plt.yticks([0.99, 1., 1.01])
+plt.xlabel("Redshift", labelpad=-0.5)
+plt.savefig("Density_evolution.png", dpi=200)

examples/Cooling/ConstantCosmoTempEvolution/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
+if [ ! -e coolingtables ]
+then
+    echo "Fetching EAGLE Cooling Tables"
+    ../getEagleCoolingTable.sh
+fi
+
+# Fetch the cooling tables
+if [ ! -e constantBox.hdf5 ]
+then
+    echo "Generating initial conditions for the uniform cosmo box example..."
+    python3 makeIC.py
+fi
+
+# Run SWIFT
+../../swift --hydro --cosmology --cooling --threads=4 const_cosmo_temp_evol.yml 2>&1 | tee output.log
+
+# Plot the result
+python3 plot_thermal_history.py cooling_box