Merge branch 'schaye08' into 'master'

Star formation following Schaye08

See merge request !705

configure.ac

@@ -1186,8 +1186,9 @@ AC_ARG_WITH([subgrid],
 with_subgrid_cooling=none
 with_subgrid_chemistry=none
 with_subgrid_tracers=none
-with_subgrid_hydro=none
+with_subgrid_entropy_floor=none
 with_subgrid_stars=none
+with_subgrid_star_formation=none
 with_subgrid_feedback=none
 
 case "$with_subgrid" in
@@ -1200,16 +1201,18 @@ case "$with_subgrid" in
 	with_subgrid_cooling=grackle
 	with_subgrid_chemistry=GEAR
 	with_subgrid_tracers=none
-	with_subgrid_hydro=gadget2
+	with_subgrid_entropy_floor=none
 	with_subgrid_stars=GEAR
+	with_subgrid_star_formation=none
 	with_subgrid_feedback=thermal
    ;;
    EAGLE)
 	with_subgrid_cooling=EAGLE
 	with_subgrid_chemistry=EAGLE
 	with_subgrid_tracers=EAGLE
-	with_subgrid_hydro=gadget2
-	with_subgrid_stars=none
+	with_subgrid_entropy_floor=EAGLE
+	with_subgrid_stars=EAGLE
+	with_subgrid_star_formation=EAGLE
 	with_subgrid_feedback=none
    ;;
    *)
@@ -1247,14 +1250,6 @@ AC_ARG_WITH([hydro],
    [with_hydro="gadget2"]
 )
 
-if test "$with_subgrid" != "none"; then
-   if test "$with_hydro" != "gadget2"; then
-      AC_MSG_ERROR([Cannot provide with-subgrid and with-hydro together])
-   else
-      with_hydro="$with_subgrid_hydro"
-   fi
-fi
-
 case "$with_hydro" in
    gadget2)
       AC_DEFINE([GADGET2_SPH], [1], [Gadget-2 SPH])
@@ -1586,7 +1581,7 @@ esac
 # Stellar model.
 AC_ARG_WITH([stars],
    [AS_HELP_STRING([--with-stars=<model>],
-      [Stellar model to use @<:@none, GEAR, debug default: none@:>@]
+      [Stellar model to use @<:@none, EAGLE, GEAR, debug default: none@:>@]
    )],
    [with_stars="$withval"],
    [with_stars="none"]
@@ -1601,10 +1596,14 @@ if test "$with_subgrid" != "none"; then
 fi
 
 case "$with_stars" in
+   EAGLE)
+      AC_DEFINE([STARS_EAGLE], [1], [EAGLE stellar model])
+   ;;
    GEAR)
       AC_DEFINE([STARS_GEAR], [1], [GEAR stellar model])
    ;;
    none)
+      AC_DEFINE([STARS_NONE], [1], [None stellar model])
    ;;
 
    *)
@@ -1682,6 +1681,62 @@ case "$with_potential" in
    ;;
 esac
 
+#  Entropy floor
+AC_ARG_WITH([entropy-floor], 
+    [AS_HELP_STRING([--with-entropy-floor=<floor>],
+       [entropy floor @<:@none, EAGLE, default: none@:>@] 
+    )],
+    [with_entropy_floor="$withval"],
+    [with_entropy_floor="none"]
+)
+if test "$with_subgrid" != "none"; then
+   if test "$with_entropy_floor" != "none"; then
+      AC_MSG_ERROR([Cannot provide with-subgrid and with-entropy-floor together])
+   else
+      with_entropy_floor="$with_subgrid_entropy_floor"
+   fi
+fi
+
+case "$with_entropy_floor" in
+   none)
+      AC_DEFINE([ENTROPY_FLOOR_NONE], [1], [No entropy floor])
+   ;;
+   EAGLE)
+      AC_DEFINE([ENTROPY_FLOOR_EAGLE], [1], [EAGLE entropy floor])
+   ;;
+   *)
+      AC_MSG_ERROR([Unknown entropy floor model])
+   ;;
+esac 
+
+#  Star formation
+AC_ARG_WITH([star-formation], 
+    [AS_HELP_STRING([--with-star-formation=<sfm>],
+       [star formation @<:@none, EAGLE, default: none@:>@] 
+    )],
+    [with_star_formation="$withval"],
+    [with_star_formation="none"]
+)
+if test "$with_subgrid" != "none"; then
+   if test "$with_star_formation" != "none"; then
+      AC_MSG_ERROR([Cannot provide with-subgrid and with-star-formation together])
+   else
+      with_star_formation="$with_subgrid_star_formation"
+   fi
+fi
+
+case "$with_star_formation" in
+   none)
+      AC_DEFINE([STAR_FORMATION_NONE], [1], [No star formation])
+   ;;
+   EAGLE)
+      AC_DEFINE([STAR_FORMATION_EAGLE], [1], [EAGLE star formation model (Schaye and Dalla Vecchia (2008))])
+   ;;
+   *)
+      AC_MSG_ERROR([Unknown star formation model])
+   ;;
+esac 
+
 #  Gravity multipole order
 AC_ARG_WITH([multipole-order],
    [AS_HELP_STRING([--with-multipole-order=<order>],
@@ -1773,11 +1828,13 @@ AC_MSG_RESULT([
    Make gravity glass  : $gravity_glass_making
    External potential  : $with_potential
 
-   Cooling function   : $with_cooling
-   Chemistry          : $with_chemistry
-   Tracers            : $with_tracers
-   Stellar model      : $with_stars
-   Feedback model     : $with_feedback
+   Entropy floor        : $with_entropy_floor
+   Cooling function     : $with_cooling
+   Chemistry            : $with_chemistry
+   Tracers              : $with_tracers
+   Stellar model        : $with_stars
+   Star formation model : $with_star_formation
+   Feedback model       : $with_feedback
 
    Individual timers           : $enable_timers
    Task debugging              : $enable_task_debugging

doc/Doxyfile.in

@@ -772,6 +772,8 @@ INPUT		       += @top_srcdir@/src/cooling/const_lambda
 INPUT		       += @top_srcdir@/src/cooling/Compton
 INPUT		       += @top_srcdir@/src/cooling/EAGLE
 INPUT		       += @top_srcdir@/src/chemistry/EAGLE
+INPUT		       += @top_srcdir@/src/entropy_floor/EAGLE
+INPUT		       += @top_srcdir@/src/star_formation/EAGLE
 INPUT		       += @top_srcdir@/src/tracers/EAGLE
 
 # This tag can be used to specify the character encoding of the source files

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

@@ -9,6 +9,66 @@ This section of the documentation gives a brief description of the
 different components of the EAGLE sub-grid model. We mostly focus on
 the parameters and values output in the snapshots.
 
+.. _EAGLE_entropy_floors:
+
+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.
+
+The first limit, labelled as ``Cool``, is typically used to prevent
+low-density high-metallicity particles to cool below the warm phase because
+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 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.
+
+.. figure:: EAGLE_entropy_floor.svg
+    :width: 400px
+    :align: center
+    :figclass: align-center
+    :alt: Phase-space diagram displaying the two entropy floors used
+	  in the EAGLE model.
+
+    Temperature-density plane with the two entropy floors used in the EAGLE
+    model indicated by the black lines. Gas particles are not allowed to be
+    below either of these two floors; they are hence forbidden to enter the
+    grey-shaded region. The floors are specified by the position in the
+    plane of the starting point of each line (black circle) and their slope
+    (dashed lines). The parameter names governing the behaviour of the
+    floors are indicated on the figure. Note that unlike what is shown on
+    the figure for clarity reasons, typical values for EAGLE runs place
+    both anchors at the same temperature.
+
+
+The model is governed by 4 parameters for each of the two
+limits. These are given in the ``EAGLEEntropyFloor`` section of the
+YAML file. The parameters are the Hydrogen number density (in
+:math:`cm^{-3}`) and temperature (in :math:`K`) of the anchor point of
+each floor as well as the power-law slope of each floor and the
+minimal over-density required to apply the limit [#f1]_. For a normal
+EAGLE run, that section of the parameter file reads:
+
+.. code:: YAML
+
+  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
+
+
 .. _EAGLE_chemical_tracers:
 
 Chemical tracers
@@ -18,7 +78,7 @@ The gas particles in the EAGLE model carry metal abundance information in the
 form of metal mass fractions. We follow explicitly 9 of the 11 elements that
 `Wiersma et al. (2009)b <http://adsabs.harvard.edu/abs/2009MNRAS.399..574W>`_
 traced in their chemical enrichment model. These are: `H`, `He`, `C`, `N`, `O`,
-`Ne`, `Mg`, `Si` and `Fe` [#f1]_. We additionally follow the total metal mass fraction
+`Ne`, `Mg`, `Si` and `Fe` [#f2]_. We additionally follow the total metal mass fraction
 (i.e. absolute metallicity) `Z`. This is typically larger than the sum of the 7
 metals that are individually traced since this will also contain the
 contribution of all the elements that are not individually followed.  We note
@@ -336,9 +396,15 @@ Black-hole accretion
 AGN feedback
 ~~~~~~~~~~~~
 
-.. [#f1] `Wiersma et al. (2009)b
+.. [#f1] Recall that in a non-cosmological run the critical density is
+	 set to 0, effectively removing the over-density
+	 constraint of the floors.
+
+.. [#f2] `Wiersma et al. (2009)b
 	 <http://adsabs.harvard.edu/abs/2009MNRAS.399..574W>`_ originally also
 	 followed explicitly `Ca` and and `S`. They are omitted in the EAGLE
 	 model but, when needed, their abundance with respect to solar is
 	 assumed to be the same as the abundance of `Si` with respect to solar
 	 (See the section :ref:`EAGLE_cooling`)
+
+

doc/RTD/source/SubgridModels/EAGLE/plot_EAGLE_entropy_floor.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,
+    "legend.fontsize": 9,
+    "xtick.labelsize": 10,
+    "ytick.labelsize": 10,
+    "text.usetex": True,
+    "figure.figsize": (3.15, 3.15),
+    "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_Jeans_rho = np.logspace(-1, 5, 1000)
+eos_Jeans_T = (eos_Jeans_rho/ 10**(-1))**(1./3.) * 4000.
+
+# 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)
+xlim(3e-8, 3e3)
+ylim(20., 2e5)
+savefig("EAGLE_entropy_floor.png", dpi=200)

examples/IsolatedGalaxy_starformation/README

+Isolated Galaxy generated by the MakeNewDisk code from Springel, Di Matteo & 
+Hernquist (2005). The done analysis in this example is similar to the work 
+done by Schaye and Dalla Vecchia (2008) (After this SD08). The default example 
+runs the simulation for a galaxy with similar mass of their fiducial model 
+and should produce plots similar to their middle pannel Figure 4. 
+
+The code can also be run for other situations to check to verify the law using
+different parameters, changes that were done in SD08 are given by:
+ - gas fraction of 10% instead of 30%, change the IC to f10.hdf5, see getIC.sh,
+   should reproduce something similar to Figure 4 left hand pannel. Requires 
+   change of fg=.1
+ - gas fraction of 90% instead of 30%, change the IC to f90.hdf5, see getIC.sh,
+   should reproduce something similar to Figure 4 right hand pannel. Requires 
+   change of fg=.9
+ - Changing the effective equation of state to adiabatic, Jeans_gamma_effective 
+   = 1.666667. Should result in something similar to Figure 5 left hand pannel
+   of SD08.
+ - Changing the effective equation of state to isothermal, Jeans_gamma_effective 
+   = 1.0000. Should result in something similar to Figure 5 middle hand pannel
+   of SD08. 
+ - Changing the slope of the Kennicutt-Schmidt law to 1.7, SchmidtLawExponent = 
+   1.7, this should result in a plot similar to Figure 6 of SD08.
+ - Increasing the density threshold by a factor of 10. thresh_norm_HpCM3 = 1.0,
+   should reproduce plot similar to Figure 7.
+ - Decreasing the density threshold by a factor of 10. thresh_norm_HpCM3 = 0.01,
+   should reproduce plot similar to Figure 7.
+ - Running with a lower resultion of a factor 8, change the IC to lowres8.hdf5,
+   see getIC.sh. 
+ - Running with a lower resultion of a factor 64, change the IC to lowres64.hdf5,
+   see getIC.sh. 
+ - Running with a lower resultion of a factor 512, change the IC to lowres512.hdf5,
+   see getIC.sh. 
+
+Other options to verify the correctness of the code is by chaning the following
+parameters:
+  - Changing the normalization to A/2 or 2A.
+  - Running the code with zero metallicity.
+  - Running the code with a factor 6 higher resolution idealized disks, use 
+    highres6.hdf5, see getIC.sh.
+  - Running with different SPH schemes like Anarchy-PU.

examples/IsolatedGalaxy_starformation/SFH.py

+"""
+Makes a movie using sphviewer and ffmpeg.
+
+Edit low_frac and up_frac to focus on a certain view of the box.
+The colour map can also be changed via colour_map.
+
+Usage: python3 makeMovie.py CoolingHalo_
+
+Written by James Willis (james.s.willis@durham.ac.uk)
+"""
+
+import glob
+import h5py as h5
+import numpy as np
+import matplotlib.pyplot as plt
+
+from tqdm import tqdm
+
+
+def getSFH(filename):
+
+    # Read the data
+    with h5.File(filename, "r") as f:
+        box_size = f["/Header"].attrs["BoxSize"][0]
+        coordinates = f["/PartType4/Coordinates"][:, :]
+        mass = f["/PartType4/Masses"][:]
+        # flag = f["/PartType4/NewStarFlag"][:]
+        birth_time = f["/PartType4/Birth_time"][:]
+
+    absmaxz = 2  # kpc
+    absmaxxy = 10  # kpc
+
+    part_mask = (
+        ((coordinates[:, 0] - box_size / 2.0) > -absmaxxy)
+        & ((coordinates[:, 0] - box_size / 2.0) < absmaxxy)
+        & ((coordinates[:, 1] - box_size / 2.0) > -absmaxxy)
+        & ((coordinates[:, 1] - box_size / 2.0) < absmaxxy)
+        & ((coordinates[:, 2] - box_size / 2.0) > -absmaxz)
+        & ((coordinates[:, 2] - box_size / 2.0) < absmaxz)
+    )  # & (flag==1)
+
+    birth_time = birth_time[part_mask]
+    mass = mass[part_mask]
+
+    histogram = np.histogram(birth_time, bins=200, weights=mass * 2e4, range=[0, 0.1])
+    values = histogram[0]
+    xvalues = (histogram[1][:-1] + histogram[1][1:]) / 2.0
+    return xvalues, values
+
+
+def getsfrsnapwide():
+
+    time = np.arange(1, 101, 1)
+    SFR_sparticles = np.zeros(100)
+    SFR_gparticles = np.zeros(100)
+    new_sparticles = np.zeros(100)
+    previous_mass = 0
+    previous_numb = 0
+    for i in tqdm(range(1, 100)):
+        # Read the data
+        filename = "output_%04d.hdf5" % i
+        with h5.File(filename, "r") as f:
+            box_size = f["/Header"].attrs["BoxSize"][0]
+            coordinates = f["/PartType0/Coordinates"][:, :]
+            SFR = f["/PartType0/SFR"][:]
+            coordinates_star = f["/PartType4/Coordinates"][:, :]
+            masses_star = f["/PartType4/Masses"][:]
+
+        absmaxz = 2  # kpc
+        absmaxxy = 10  # kpc
+
+        part_mask = (
+            ((coordinates[:, 0] - box_size / 2.0) > -absmaxxy)
+            & ((coordinates[:, 0] - box_size / 2.0) < absmaxxy)
+            & ((coordinates[:, 1] - box_size / 2.0) > -absmaxxy)
+            & ((coordinates[:, 1] - box_size / 2.0) < absmaxxy)
+            & ((coordinates[:, 2] - box_size / 2.0) > -absmaxz)
+            & ((coordinates[:, 2] - box_size / 2.0) < absmaxz)
+            & (SFR > 0)
+        )
+
+        SFR = SFR[part_mask]
+
+        total_SFR = np.sum(SFR)
+        SFR_gparticles[i] = total_SFR * 10
+
+    return time[:-1], SFR_gparticles[1:]
+
+
+if __name__ == "__main__":
+
+    time, SFR1 = getsfrsnapwide()  # , SFR2, SFR_error = getsfrsnapwide()
+    time2, SFR3 = getSFH("output_%04d.hdf5" % 100)
+    plt.plot(time2[1:] * 1e3, SFR3[1:], label="Using birth_time of star particles")
+    plt.plot(time, SFR1, label="Using SFR of gas particles", color="g")
+    plt.xlabel("Time (Myr)")
+    plt.ylabel("SFH ($\\rm M_\odot \\rm yr^{-1}$)")
+    plt.ylim(0, 20)
+    plt.legend()
+    plt.savefig("SFH.png")

examples/IsolatedGalaxy_starformation/getIC.sh

+#!/bin/bash 
+
+wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/fid.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/f10.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/f90.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/lowres8.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/lowres64.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/lowres512.hdf5
+# wget https://www.strw.leidenuniv.nl/~nobels/swiftdata/highres6.hdf5

examples/IsolatedGalaxy_starformation/isolated_galaxy.yml

+# Define the system of units to use internally.
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.9891E43   # 10^10 solar masses 
+  UnitLength_in_cgs:   3.08567758E21   # 1 kpc 
+  UnitVelocity_in_cgs: 1E5   # km/s
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters for the self-gravity scheme
+Gravity:
+  mesh_side_length:       32        # Number of cells along each axis for the periodic gravity mesh.
+  eta:          0.025               # Constant dimensionless multiplier for time integration.
+  theta:        0.7                 # Opening angle (Multipole acceptance criterion).
+  comoving_softening:     0.01 # Comoving softening length (in internal units).
+  max_physical_softening: 0.01    # Physical softening length (in internal units).
+
+# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
+TimeIntegration:
+  time_begin:        0.    # The starting time of the simulation (in internal units).
+  time_end:          0.1    # The end time 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
+Snapshots:
+  basename:   output      # Common part of the name of output files
+  time_first: 0.          # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
+  delta_time: 0.001        # Time difference between consecutive outputs (in internal units)
+
+  
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:           1e-2     # Time between statistics output
+  time_first:             0.     # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units)
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  fid.hdf5 # The file to read
+  periodic:                    0    # Are we running with periodic ICs?
+
+# 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.
+  h_tolerance:           1e-4     # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
+  h_max:                 10.      # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
+  max_volume_change:     1.4      # (Optional) Maximal allowed change of kernel volume over one time-step.
+  max_ghost_iterations:  30       # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
+  minimal_temperature:   100        # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
+  H_ionization_temperature: 1e4   # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas.
+
+EAGLECooling:
+  dir_name:                ./coolingtables/  # Location of the Wiersma+08 cooling tables
+  H_reion_z:               11.5              # Redshift of Hydrogen re-ionization
+  He_reion_z_centre:       3.5               # Redshift of the centre of the Helium re-ionization Gaussian
+  He_reion_z_sigma:        0.5               # Spread in redshift of the  Helium re-ionization Gaussian
+  He_reion_eV_p_H:         2.0               # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
+
+# Primordial abundances
+EAGLEChemistry:
+  init_abundance_metal:     0.0129          # Inital fraction of particle mass in *all* metals
+  init_abundance_Hydrogen:  0.7065       # Inital fraction of particle mass in Hydrogen
+  init_abundance_Helium:    0.2806        # Inital fraction of particle mass in Helium
+  init_abundance_Carbon:    0.00207        # Inital fraction of particle mass in Carbon
+  init_abundance_Nitrogen:  0.000836        # Inital fraction of particle mass in Nitrogen
+  init_abundance_Oxygen:    0.00549        # Inital fraction of particle mass in Oxygen
+  init_abundance_Neon:      0.00141        # Inital fraction of particle mass in Neon
+  init_abundance_Magnesium: 0.000591        # Inital fraction of particle mass in Magnesium
+  init_abundance_Silicon:   0.000683        # Inital fraction of particle mass in Silicon
+  init_abundance_Iron:      0.0011        # Inital fraction of particle mass in Iron
+
+# Hernquist potential parameters