diff --git a/doc/RTD/source/Cooling/index.rst b/doc/RTD/source/Cooling/index.rst
deleted file mode 100644
index 00b84489c647bbca4948714a95b3d5cd5fb3cce1..0000000000000000000000000000000000000000
--- a/doc/RTD/source/Cooling/index.rst
+++ /dev/null
@@ -1,80 +0,0 @@
-.. Equation of State
- Loic Hausammann, 7th April 2018
-
-.. _cooling:
-
-Cooling
-=======
-
-Currently, we have 5 different cooling (EAGLE, Grackle, const-lambda, const-du
-and none). Three of them are easily solved analytically (const-lambda,
-const-du and none) while the two last requires complex chemical networks.
-
-
-Equations
----------
-
-The first table compares the different analytical cooling while the next ones
-are specific to a given cooling. The quantities are the internal energy (\\( u
-\\)), the density \\( rho \\), the element mass fraction (\\( X_i \\)), the
-cooling function (\\(\\Lambda\\), the proton mass (\\( m_H \\)) and the time
-step condition (\\( t\_\\text{step}\\)). If not specified otherwise, all
-cooling contains a temperature floor avoiding negative temperature.
-
-.. csv-table:: Analytical Cooling
- :header: "Variable", "Const-Lambda", "Const-du", "None"
-
- "\\( \\frac{ \\mathrm{d}u }{ \\mathrm{d}t } \\)", "\\( -\\Lambda \\frac{\\rho^2 X_H^2}{\\rho m_H^2} \\)", "const", "0"
- "\\( \\Delta t\_\\text{max} \\)", "\\( t\_\\text{step} \\frac{u}{\\left|\\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \\)", "\\( t\_\\text{step} \\frac{u}{\\ \\left| \\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \\)", "None"
-
-
-Grackle
-~~~~~~~
-
-Grackle is a chemistry and cooling library presented in `B. Smith et al. 2016 <https://arxiv.org/abs/1610.09591>`_
-(do not forget to cite if used). Four different modes are available:
-equilibrium, 6 species network (H, H\\( ^+ \\), e\\( ^- \\), He, He\\( ^+ \\)
-and He\\( ^{++} \\)), 9 species network (adds H\\(^-\\), H\\(_2\\) and
-H\\(_2^+\\)) and 12 species (adds D, D\\(^+\\) and HD). Following the same
-order, the swift cooling options are ``grackle``, ``grackle1``, ``grackle2``
-and ``grackle3`` (the numbers correspond to the value of
-``primordial_chemistry`` in Grackle). It also includes some self-shielding
-methods and UV background. In order to use the Grackle cooling, you will need
-to provide an HDF5 table computed by Cloudy.
-
-When starting a simulation without providing the different fractions, the code
-supposes an equilibrium and computes the fractions automatically.
-
-In order to compile SWIFT with Grackle, you need to provide the options ``with-grackle``
-and ``with-chemistry``.
-
-You will need a Grackle version later than 3.1. To compile it, run
-the following commands from the root directory of Grackle:
-``./configure; cd src/clib``.
-Update the variables ``LOCAL_HDF5_INSTALL`` and ``MACH_INSTALL_PREFIX`` in
-the file ``src/clib/Make.mach.linux-gnu``.
-Finish with ``make machine-linux-gnu; make && make install``.
-If you encounter any problem, you can look at the `Grackle documentation <https://grackle.readthedocs.io/en/latest/>`_
-
-You can now provide the path given for ``MACH_INSTALL_PREFIX`` to ``with-grackle``.
-
-Eagle
-~~~~~
-
-TODO
-
-How to Implement a New Cooling
-------------------------------
-
-The developer should provide at least one function for:
- * writing the cooling name in HDF5
- * cooling a particle
- * the maximal time step possible
- * initializing a particle
- * computing the total energy radiated by a particle
- * initializing the cooling parameters
- * printing the cooling type
-
-For implementation details, see ``src/cooling/none/cooling.h``
-
-See :ref:`new_option` for the full list of changes required.
diff --git a/doc/RTD/source/SubgridModels/Basic/index.rst b/doc/RTD/source/SubgridModels/Basic/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..031a1ed61bb1cae4f916b58b926f5269b10c7057
--- /dev/null
+++ b/doc/RTD/source/SubgridModels/Basic/index.rst
@@ -0,0 +1,51 @@
+.. Basic sub-grid model
+ Matthieu Schaller, 20th December 2018
+
+
+Basic model
+===========
+
+
+Cooling: Analytic models
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+Currently, we have 3 different simple cooling models (const-lambda, const-du
+and Compton). These are all based on analytic formulas and can be used
+to quickly understand how the cooling interacts with the rest of the
+code before moving to more complex models.
+
+Equations
+---------
+
+The first table compares the different analytical cooling while the next ones
+are specific to a given cooling. The quantities are the internal energy (\\( u
+\\)), the density \\( rho \\), the element mass fraction (\\( X_i \\)), the
+cooling function (\\(\\Lambda\\), the proton mass (\\( m_H \\)) and the time
+step condition (\\( t\_\\text{step}\\)). If not specified otherwise, all
+cooling contains a temperature floor avoiding negative temperature.
+
+.. csv-table:: Analytical Cooling
+ :header: "Variable", "Const-Lambda", "Const-du"
+
+ "\\( \\frac{ \\mathrm{d}u }{ \\mathrm{d}t } \\)", "\\( -\\Lambda \\frac{\\rho^2 X_H^2}{\\rho m_H^2} \\)", "const"
+ "\\( \\Delta t\_\\text{max} \\)", "\\( t\_\\text{step} \\frac{u}{\\left|\\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \\)", "\\( t\_\\text{step} \\frac{u}{\\ \\left| \\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \\)"
+
+TODO: Add description of the parameters and units.
+
+TODO: Add Compton cooling model
+
+How to Implement a New Cooling
+------------------------------
+
+The developer should provide at least one function for:
+ * writing the cooling name in HDF5
+ * cooling a particle
+ * the maximal time step possible
+ * initializing a particle
+ * computing the total energy radiated by a particle
+ * initializing the cooling parameters
+ * printing the cooling type
+
+For implementation details, see ``src/cooling/none/cooling.h``
+
+See :ref:`new_option` for the full list of changes required.
diff --git a/doc/RTD/source/SubgridModels/EAGLE/index.rst b/doc/RTD/source/SubgridModels/EAGLE/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..332524240b3d8933d6425c369ed1fc5ec9247a1b
--- /dev/null
+++ b/doc/RTD/source/SubgridModels/EAGLE/index.rst
@@ -0,0 +1,80 @@
+.. EAGLE sub-grid model
+ Matthieu Schaller, 20th December 2018
+
+
+EAGLE model
+===========
+
+Chemical tracers
+~~~~~~~~~~~~~~~~
+
+The gas particles in the EAGLE model carry metal abundance information
+in the form of metal mass fractions. We follow the following 9
+elements: `H`, `He`, `C`, `N`, `O`, `Ne`, `Mg`, `Si` and `Fe`. We
+additionally follow the total metal mass fraction `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.
+
+As part of the diagnostics, we additionally trace the elements coming
+from the different stellar evolution channels. We store for each
+particle the total mass coming from all the SNIa that enriched that
+particle and the metal mass fraction from SNIa. This is the fraction
+of the *total* gas mass that is in the form of metals originating from
+SNIa stars. By construction this fraction will be smaller than the
+total metal mass fraction. The same tracers exist for the SNII and AGB
+channels. Finally, we also compute the iron gas fraction from
+SNIa. This it the fraction of the *total* gas mass that is made of
+iron originating from SNIa explosions.
+
+We finally also compute the smoothed versions of the individual
+element mass fractions, of the total metal mass fractions, and of the
+iron gas fraction from SNIa.
+
+The chemistry module in ``src/chemistry/EAGLE`` includes all the arrays
+that are added to the particles and the functions used to compute the
+smoothed elements.
+
+In the snapshots, we output:
+
++------------------------------+--------------------------------+-----------+-----------------------------+
+| Name | Description | Units | Comments |
++==============================+================================+===========+=============================+
+| ``ElementAbundance`` | | Fraction of the gas mass | [-] | | Array of length |
+| | | in the different elements | | | 9 for each particle |
++------------------------------+--------------------------------+-----------+-----------------------------+
+| ``SmoothedElementAbundance`` | | Fraction of the gas mass | [-] | | Array of length |
+| | | in the different elements | | | 9 for each particle |
+| | | smoothed over SPH neighbours | | |
++------------------------------+--------------------------------+-----------+-----------------------------+
+| ``Metallicity`` | | Fraction of the gas mass | [-] | |
+| | | in *all* metals | | |
++------------------------------+--------------------------------+-----------+-----------------------------+
+| ``SmoothedMetallicity`` | | Fraction of the gas mass | [-] | |
+| | | in *all* metals | | |
+| | | smoothed over SPH neighbours | | |
++------------------------------+--------------------------------+-----------+-----------------------------+
+
+Cooling: Wiersma+2008a
+~~~~~~~~~~~~~~~~~~~~~~
+
+Particle tracers
+~~~~~~~~~~~~~~~~
+
+Star formation: Schaye+2008
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Stellar enrichment: Wiersma+2008b
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Supernova feedback: Schaye+2012
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Black-hole creation
+~~~~~~~~~~~~~~~~~~~
+
+Black-hole accretion
+~~~~~~~~~~~~~~~~~~~~
+
+AGN feedback
+~~~~~~~~~~~~
diff --git a/doc/RTD/source/SubgridModels/GEAR/index.rst b/doc/RTD/source/SubgridModels/GEAR/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..2a211759bfc4895fd07279b72f78200d6ea47546
--- /dev/null
+++ b/doc/RTD/source/SubgridModels/GEAR/index.rst
@@ -0,0 +1,37 @@
+.. GEAR sub-grid model
+ Matthieu Schaller, 20th December 2018
+
+
+GEAR model
+===========
+
+
+Cooling: Grackle
+~~~~~~~~~~~~~~~~
+
+Grackle is a chemistry and cooling library presented in `B. Smith et al. 2016 <https://arxiv.org/abs/1610.09591>`_
+(do not forget to cite if used). Four different modes are available:
+equilibrium, 6 species network (H, H\\( ^+ \\), e\\( ^- \\), He, He\\( ^+ \\)
+and He\\( ^{++} \\)), 9 species network (adds H\\(^-\\), H\\(_2\\) and
+H\\(_2^+\\)) and 12 species (adds D, D\\(^+\\) and HD). Following the same
+order, the swift cooling options are ``grackle``, ``grackle1``, ``grackle2``
+and ``grackle3`` (the numbers correspond to the value of
+``primordial_chemistry`` in Grackle). It also includes some self-shielding
+methods and UV background. In order to use the Grackle cooling, you will need
+to provide an HDF5 table computed by Cloudy.
+
+When starting a simulation without providing the different fractions, the code
+supposes an equilibrium and computes the fractions automatically.
+
+In order to compile SWIFT with Grackle, you need to provide the options ``with-grackle``
+and ``with-chemistry``.
+
+You will need a Grackle version later than 3.1. To compile it, run
+the following commands from the root directory of Grackle:
+``./configure; cd src/clib``.
+Update the variables ``LOCAL_HDF5_INSTALL`` and ``MACH_INSTALL_PREFIX`` in
+the file ``src/clib/Make.mach.linux-gnu``.
+Finish with ``make machine-linux-gnu; make && make install``.
+If you encounter any problem, you can look at the `Grackle documentation <https://grackle.readthedocs.io/en/latest/>`_
+
+You can now provide the path given for ``MACH_INSTALL_PREFIX`` to ``with-grackle``.
diff --git a/doc/RTD/source/SubgridModels/index.rst b/doc/RTD/source/SubgridModels/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..c1b8bc858c527f0ec6834e8de163306dd1be66cc
--- /dev/null
+++ b/doc/RTD/source/SubgridModels/index.rst
@@ -0,0 +1,18 @@
+.. Subgrid Models
+ Matthieu Schaller, 20th December 2018
+
+.. _subgrid:
+
+Galaxy Formation Subgrid Models
+===============================
+
+Multiple models are available in SWIFT. The 'Basic' model can
+be use as an empty canvas to be copied to create additional models.
+
+.. toctree::
+ :maxdepth: 2
+ :caption: Available models:
+
+ Basic/index
+ EAGLE/index
+ GEAR/index
diff --git a/doc/RTD/source/index.rst b/doc/RTD/source/index.rst
index d148398c1bd77eafbce5e0037457b34efddb4eca..b9370c3f24b2ffb3c5174f2fe99fb9ec610e18f6 100644
--- a/doc/RTD/source/index.rst
+++ b/doc/RTD/source/index.rst
@@ -19,7 +19,7 @@ difference is the parameter file that will need to be adapted for SWIFT.
ParameterFiles/index
InitialConditions/index
HydroSchemes/index
- Cooling/index
+ SubgridModels/index
EquationOfState/index
ExternalPotentials/index
NewOption/index
diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml
index df2c29c0b612eff377423b7bb76e2c8e1e530df1..6207eee94429411b08cc55c6484b7e679aec1092 100644
--- a/examples/CoolingBox/coolingBox.yml
+++ b/examples/CoolingBox/coolingBox.yml
@@ -70,8 +70,6 @@ EAGLEChemistry:
InitAbundance_Magnesium: 5.907064e-4
InitAbundance_Silicon: 6.825874e-4
InitAbundance_Iron: 1.1032152e-3
- CalciumOverSilicon: 0.0941736
- SulphurOverSilicon: 0.6054160
GearChemistry:
InitialMetallicity: 0.01295
diff --git a/examples/CoolingRates/cooling_rates.yml b/examples/CoolingRates/cooling_rates.yml
index e0ac9f691cf64b292c50a36d2b1878bf3a368975..2bda9b25c40046a5363ff8ebcc08aa82ae10b798 100644
--- a/examples/CoolingRates/cooling_rates.yml
+++ b/examples/CoolingRates/cooling_rates.yml
@@ -26,9 +26,6 @@ EAGLEChemistry:
InitAbundance_Magnesium: 5.907064e-4
InitAbundance_Silicon: 6.825874e-4
InitAbundance_Iron: 1.1032152e-3
- CalciumOverSilicon: 0.0941736
- SulphurOverSilicon: 0.6054160
-
EagleCooling:
filename: ./coolingtables/
diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index 1d19320491130ce6d98ac8cd0adff8753f2c2e54..430fc6610881984ecf5ec1b4e3fc6f92ee065eae 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -68,8 +68,6 @@ EAGLEChemistry: # Solar abundances
InitAbundance_Magnesium: 5.907064e-4
InitAbundance_Silicon: 6.825874e-4
InitAbundance_Iron: 1.1032152e-3
- CalciumOverSilicon: 0.0941736
- SulphurOverSilicon: 0.6054160
EagleCooling:
filename: ./coolingtables/
diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml
index 3e1e5907df7bf76e8ad96a1c9d3c667945c090c8..11401699c308bdf22832708f4a208edc2ca572a8 100644
--- a/examples/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_25/eagle_25.yml
@@ -77,8 +77,6 @@ EAGLEChemistry: # Solar abundances
InitAbundance_Magnesium: 5.907064e-4
InitAbundance_Silicon: 6.825874e-4
InitAbundance_Iron: 1.1032152e-3
- CalciumOverSilicon: 0.0941736
- SulphurOverSilicon: 0.6054160
EagleCooling:
filename: ./coolingtables/
diff --git a/examples/EAGLE_50/eagle_50.yml b/examples/EAGLE_50/eagle_50.yml
index 87bef197cd357116d4b67015c5f58774c23d36b7..91c9ece9b2c2d3ff01ba831efa4fbb5a954a9537 100644
--- a/examples/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_50/eagle_50.yml
@@ -70,8 +70,6 @@ EAGLEChemistry: # Solar abundances
InitAbundance_Magnesium: 5.907064e-4
InitAbundance_Silicon: 6.825874e-4
InitAbundance_Iron: 1.1032152e-3
- CalciumOverSilicon: 0.0941736
- SulphurOverSilicon: 0.6054160
EagleCooling:
filename: ./coolingtables/
diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml
index 95a5d3398f13d42700963a0abbae118083662440..3e61666bad2e982bb01cabae5a8af9952a027772 100644
--- a/examples/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_6/eagle_6.yml
@@ -81,8 +81,6 @@ EAGLEChemistry: # Solar abundances
InitAbundance_Magnesium: 5.907064e-4
InitAbundance_Silicon: 6.825874e-4
InitAbundance_Iron: 1.1032152e-3
- CalciumOverSilicon: 0.0941736
- SulphurOverSilicon: 0.6054160
EagleCooling:
filename: ./coolingtables/
diff --git a/examples/SantaBarbara/santa_barbara.yml b/examples/SantaBarbara/santa_barbara.yml
index f217fc6a67db92eb9ea5e4d541343e7bf89c5860..964d0d448c94803857ca3e4e0d6ab43ad71488c8 100644
--- a/examples/SantaBarbara/santa_barbara.yml
+++ b/examples/SantaBarbara/santa_barbara.yml
@@ -70,8 +70,6 @@ EAGLEChemistry:
InitAbundance_Magnesium: 0.0
InitAbundance_Silicon: 0.0
InitAbundance_Iron: 0.0
- CalciumOverSilicon: 0.0
- SulphurOverSilicon: 0.0
EagleCooling:
filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/
diff --git a/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml b/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml
index 969626d97bf98334dd87cb9dc6862b77795b1643..0fdad57cb194476cbbea99d4db3d80591c66cae9 100644
--- a/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml
@@ -82,5 +82,3 @@ EAGLEChemistry:
InitAbundance_Magnesium: 0.
InitAbundance_Silicon: 0.
InitAbundance_Iron: 0.
- CalciumOverSilicon: 0.
- SulphurOverSilicon: 0.
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 6258782ab802ae85d783543bdeedf34f538333cf..11cd2854e39683cee4859be16fc6964327cab240 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -248,12 +248,15 @@ LambdaCooling:
lambda_nH2_cgs: 1e-22 # 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.
+# Parameters of the EAGLE cooling model (Wiersma+08 cooling tables).
EagleCooling:
- filename: ./coolingtables/ # Location of the Wiersma+08 cooling tables
- reionisation_redshift: 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_pH: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
+ filename: ./coolingtables/ # Location of the Wiersma+08 cooling tables
+ reionisation_redshift: 8.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_pH: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
+ CalciumOverSiliconInSolar: 1. # (Optional) Ratio of Ca/Si to use in units of solar. If set to 1, the code uses [Ca/Si] = 0, i.e. Ca/Si = 0.0941736.
+ SulphurOverSiliconInSolar: 1. # (Optional) Ratio of S/Si to use in units of solar. If set to 1, the code uses [S/Si] = 0, i.e. S/Si = 0.6054160.
# Cooling with Grackle 3.0
GrackleCooling:
@@ -282,8 +285,6 @@ EAGLEChemistry:
InitAbundance_Magnesium: 0.000 # Inital fraction of particle mass in Magnesium
InitAbundance_Silicon: 0.000 # Inital fraction of particle mass in Silicon
InitAbundance_Iron: 0.000 # Inital fraction of particle mass in Iron
- CalciumOverSilicon: 0.0941736 # Constant ratio of Calcium over Silicon abundance
- SulphurOverSilicon: 0.6054160 # Constant ratio of Sulphur over Silicon abundance
# Structure finding options (requires velociraptor)
StructureFinding:
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index 41a62c2d9f1387e58fe9027c4ca7ce0dee144514..9e9dc3b3232b326c5d6dc9487a75dce3b82dd401 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -188,12 +188,6 @@ static INLINE void chemistry_init_backend(struct swift_params* parameter_file,
data->initial_metal_mass_fraction[elem] =
parser_get_param_float(parameter_file, buffer);
}
-
- /* Read the constant ratios */
- data->calcium_over_silicon_ratio = parser_get_param_float(
- parameter_file, "EAGLEChemistry:CalciumOverSilicon");
- data->sulphur_over_silicon_ratio = parser_get_param_float(
- parameter_file, "EAGLEChemistry:SulphurOverSilicon");
}
}
diff --git a/src/chemistry/EAGLE/chemistry_struct.h b/src/chemistry/EAGLE/chemistry_struct.h
index 9093709e62d0af638ae485bd1154a4537791e84a..f5e47347f8b6a910640624ddfc0b5968242eedf6 100644
--- a/src/chemistry/EAGLE/chemistry_struct.h
+++ b/src/chemistry/EAGLE/chemistry_struct.h
@@ -45,12 +45,6 @@ struct chemistry_global_data {
/*! Fraction of the particle mass in *all* metals at the start of the run */
float initial_metal_mass_fraction_total;
-
- /*! Constant ratio of Calcium over Silicium */
- float calcium_over_silicon_ratio;
-
- /*! Constant ratio of Calcium over Silicium */
- float sulphur_over_silicon_ratio;
};
/**
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index b94d60ff4346d52efd37df0795b21133e76c925b..31ab022fc934c77db6ac0371396acc3204cd0f56 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -466,7 +466,10 @@ void cooling_cool_part(const struct phys_const *restrict phys_const,
(See cosmology theory document for the derivation) */
const double delta_redshift = -dt * cosmo->H * cosmo->a_inv;
- /* Get this particle's abundance ratios */
+ /* Get this particle's abundance ratios compared to solar
+ * Note that we need to add S and Ca that are in the tables but not tracked
+ * by the particles themselves.
+ * The order is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe] */
float abundance_ratio[chemistry_element_count + 2];
abundance_ratio_to_solar(p, cooling, abundance_ratio);
@@ -739,10 +742,10 @@ void cooling_init_backend(struct swift_params *parameter_file,
cooling->cooling_table_path);
cooling->H_reion_z = parser_get_param_float(
parameter_file, "EagleCooling:reionisation_redshift");
- cooling->calcium_over_silicon_ratio = parser_get_param_float(
- parameter_file, "EAGLEChemistry:CalciumOverSilicon");
- cooling->sulphur_over_silicon_ratio = parser_get_param_float(
- parameter_file, "EAGLEChemistry:SulphurOverSilicon");
+ cooling->Ca_over_Si_ratio_in_solar = parser_get_opt_param_float(
+ parameter_file, "EAGLECooling:CalciumOverSiliconInSolar", 1.f);
+ cooling->S_over_Si_ratio_in_solar = parser_get_opt_param_float(
+ parameter_file, "EAGLECooling:SulphurOverSiliconInSolar", 1.f);
cooling->He_reion_z_centre =
parser_get_param_float(parameter_file, "EagleCooling:He_reion_z_centre");
cooling->He_reion_z_sigma =
@@ -866,6 +869,7 @@ void cooling_clean(struct cooling_function_data *cooling) {
free(cooling->HeFrac);
free(cooling->Therm);
free(cooling->SolarAbundances);
+ free(cooling->SolarAbundances_inv);
/* Free the tables */
free(cooling->table.metal_heating);
diff --git a/src/cooling/EAGLE/cooling_rates.h b/src/cooling/EAGLE/cooling_rates.h
index c78be6a30849a7b8a79e56d7fbdefe2eede866f3..eb9e985d6588eaa69b15b0626e60a6a160db5cd7 100644
--- a/src/cooling/EAGLE/cooling_rates.h
+++ b/src/cooling/EAGLE/cooling_rates.h
@@ -29,50 +29,68 @@
#include "interpolate.h"
/**
- * @brief Calculate ratio of particle element abundances
- * to solar abundance.
-
- * Multiple if statements are necessary because order of elements
- * in tables is different from chemistry_element enum.
- * Tables: H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe
- * Enum: H, He, C, N, O, Ne, Mg, Si, Fe
- * The order in ratio_solar is:
- * H, He, C, N, O, Ne, Mg, Si, Fe, S, Ca
- * Hence Fe, S, Ca need to be treated separately to be put in the
- * correct place in the output array.
- *
- * @param p Pointer to #part struct
- * @param cooling #cooling_function_data struct
- * @param ratio_solar Pointer to array or ratios
+ * @brief Compute ratio of mass fraction to solar mass fraction
+ * for each element carried by a given particle.
+ *
+ * The solar abundances are taken from the tables themselves.
+ *
+ * The EAGLE chemistry model does not track S and Ca. We assume
+ * that their abundance with respect to solar is the same as
+ * the ratio for Si.
+ * We optionally apply a correction if the user asked for a different
+ * ratio.
+ *
+ * We also re-order the elements such that they match the order of the
+ * tables. This is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe].
+ *
+ * The solar abundances table is arranged as
+ * [H, He, C, N, O, Ne, Mg, Si, S, C, Fe].
+ *
+ * @param p Pointer to #part struct.
+ * @param cooling #cooling_function_data struct.
+ * @param ratio_solar (return) Array of ratios to solar abundances.
*/
__attribute__((always_inline)) INLINE void abundance_ratio_to_solar(
const struct part *p, const struct cooling_function_data *cooling,
float ratio_solar[chemistry_element_count + 2]) {
- /* compute ratios for all elements */
- for (enum chemistry_element elem = chemistry_element_H;
- elem < chemistry_element_count; elem++) {
- if (elem == chemistry_element_Fe) {
- /* NOTE: solar abundances have iron last with calcium and sulphur directly
- * before, hence +2 */
- ratio_solar[elem] = p->chemistry_data.metal_mass_fraction[elem] /
- cooling->SolarAbundances[elem + 2];
- } else {
- ratio_solar[elem] = p->chemistry_data.metal_mass_fraction[elem] /
- cooling->SolarAbundances[elem];
- }
- }
+ ratio_solar[0] = p->chemistry_data.metal_mass_fraction[chemistry_element_H] *
+ cooling->SolarAbundances_inv[0];
+
+ ratio_solar[1] = p->chemistry_data.metal_mass_fraction[chemistry_element_He] *
+ cooling->SolarAbundances_inv[1];
+
+ ratio_solar[2] = p->chemistry_data.metal_mass_fraction[chemistry_element_C] *
+ cooling->SolarAbundances_inv[2];
+
+ ratio_solar[3] = p->chemistry_data.metal_mass_fraction[chemistry_element_N] *
+ cooling->SolarAbundances_inv[3];
- /* assign ratios for Ca and S, note positions of these elements occur before
- * Fe */
- ratio_solar[chemistry_element_count] =
- p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
- cooling->sulphur_over_silicon_ratio /
- cooling->SolarAbundances[chemistry_element_count - 1];
- ratio_solar[chemistry_element_count + 1] =
- p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
- cooling->calcium_over_silicon_ratio /
- cooling->SolarAbundances[chemistry_element_count];
+ ratio_solar[4] = p->chemistry_data.metal_mass_fraction[chemistry_element_O] *
+ cooling->SolarAbundances_inv[4];
+
+ ratio_solar[5] = p->chemistry_data.metal_mass_fraction[chemistry_element_Ne] *
+ cooling->SolarAbundances_inv[5];
+
+ ratio_solar[6] = p->chemistry_data.metal_mass_fraction[chemistry_element_Mg] *
+ cooling->SolarAbundances_inv[6];
+
+ ratio_solar[7] = p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
+ cooling->SolarAbundances_inv[7];
+
+ /* For S, we use the same ratio as Si */
+ ratio_solar[8] = p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
+ cooling->SolarAbundances_inv[7] *
+ cooling->S_over_Si_ratio_in_solar;
+
+ /* For Ca, we use the same ratio as Si */
+ ratio_solar[9] = p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
+ cooling->SolarAbundances_inv[7] *
+ cooling->Ca_over_Si_ratio_in_solar;
+
+ ratio_solar[10] =
+ p->chemistry_data.metal_mass_fraction[chemistry_element_Fe] *
+ cooling->SolarAbundances_inv[10];
}
/**
@@ -585,7 +603,7 @@ INLINE static double eagle_metal_cooling_rate(
#endif
}
- const double abundance_ratio =
+ const double electron_abundance_ratio =
H_plus_He_electron_abundance / solar_electron_abundance;
/**********************/
@@ -596,22 +614,28 @@ INLINE static double eagle_metal_cooling_rate(
* electron abundance to solar electron abundance then by the ratio of the
* particle metal abundance to solar metal abundance. */
- double lambda_metal[eagle_cooling_N_metal];
+ double lambda_metal[eagle_cooling_N_metal + 2] = {0.};
if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) {
- for (int elem = 0; elem < eagle_cooling_N_metal; elem++) {
+ /* Loop over the metals (ignore H and He) */
+ for (int elem = 2; elem < eagle_cooling_N_metal + 2; elem++) {
+
+ if (solar_ratio[elem] > 0.) {
- lambda_metal[elem] =
- interpolation_3d_no_x(cooling->table.metal_heating, /* */
- elem, n_H_index, T_index, /* */
- /*delta_elem=*/0.f, d_n_H, d_T, /* */
- eagle_cooling_N_metal, /* */
- eagle_cooling_N_density, /* */
- eagle_cooling_N_temperature); /* */
+ /* Note that we do not interpolate along the x-axis
+ * (element dimension) */
+ lambda_metal[elem] =
+ interpolation_3d_no_x(cooling->table.metal_heating, /* */
+ elem - 2, n_H_index, T_index, /* */
+ /*delta_elem=*/0.f, d_n_H, d_T, /* */
+ eagle_cooling_N_metal, /* */
+ eagle_cooling_N_density, /* */
+ eagle_cooling_N_temperature); /* */
- lambda_metal[elem] *= abundance_ratio;
- lambda_metal[elem] *= solar_ratio[elem + 2];
+ lambda_metal[elem] *= electron_abundance_ratio;
+ lambda_metal[elem] *= solar_ratio[elem];
+ }
#ifdef TO_BE_DONE
/* compute values at temperature gridpoints above and below input
@@ -636,19 +660,25 @@ INLINE static double eagle_metal_cooling_rate(
} else {
- for (int elem = 0; elem < eagle_cooling_N_metal; elem++) {
+ /* Loop over the metals (ignore H and He) */
+ for (int elem = 2; elem < eagle_cooling_N_metal + 2; elem++) {
- lambda_metal[elem] = interpolation_4d_no_x(
- cooling->table.metal_heating, /* */
- elem, /*z_index=*/0, n_H_index, T_index, /* */
- /*delta_elem=*/0.f, cooling->dz, d_n_H, d_T, /* */
- eagle_cooling_N_metal, /* */
- eagle_cooling_N_loaded_redshifts, /* */
- eagle_cooling_N_density, /* */
- eagle_cooling_N_temperature); /* */
+ if (solar_ratio[elem] > 0.) {
- lambda_metal[elem] *= abundance_ratio;
- lambda_metal[elem] *= solar_ratio[elem + 2];
+ /* Note that we do not interpolate along the x-axis
+ * (element dimension) */
+ lambda_metal[elem] = interpolation_4d_no_x(
+ cooling->table.metal_heating, /* */
+ elem - 2, /*z_index=*/0, n_H_index, T_index, /* */
+ /*delta_elem=*/0.f, cooling->dz, d_n_H, d_T, /* */
+ eagle_cooling_N_metal, /* */
+ eagle_cooling_N_loaded_redshifts, /* */
+ eagle_cooling_N_density, /* */
+ eagle_cooling_N_temperature); /* */
+
+ lambda_metal[elem] *= electron_abundance_ratio;
+ lambda_metal[elem] *= solar_ratio[elem];
+ }
#ifdef TO_BE_DONE
/* compute values at temperature gridpoints above and below input
@@ -675,14 +705,14 @@ INLINE static double eagle_metal_cooling_rate(
}
if (element_lambda != NULL) {
- for (int elem = 0; elem < eagle_cooling_N_metal; ++elem) {
- element_lambda[elem + 2] = lambda_metal[elem];
+ for (int elem = 2; elem < eagle_cooling_N_metal + 2; ++elem) {
+ element_lambda[elem] = lambda_metal[elem];
}
}
/* Sum up all the contributions */
double Lambda_net = Lambda_free + Lambda_Compton;
- for (int elem = 0; elem < eagle_cooling_N_metal; ++elem) {
+ for (int elem = 2; elem < eagle_cooling_N_metal + 2; ++elem) {
Lambda_net += lambda_metal[elem];
}
diff --git a/src/cooling/EAGLE/cooling_struct.h b/src/cooling/EAGLE/cooling_struct.h
index 09a16e8c10e5e0a1f03ea3455e3ffff8c942a982..ae18d8b864a3e41529831169bf0d15b0e917b50b 100644
--- a/src/cooling/EAGLE/cooling_struct.h
+++ b/src/cooling/EAGLE/cooling_struct.h
@@ -65,20 +65,23 @@ struct cooling_function_data {
/*! Internal energy bins */
float *Therm;
- /*! Solar mass fractions */
+ /*! Mass fractions of elements for solar abundances (from the tables) */
float *SolarAbundances;
+ /*! Inverse of the solar mass fractions */
+ float *SolarAbundances_inv;
+
/*! Filepath to the directory containing the HDF5 cooling tables */
char cooling_table_path[eagle_table_path_name_length];
/*! Redshit of H reionization */
float H_reion_z;
- /*! Ca over Si abundance ratio */
- float calcium_over_silicon_ratio;
+ /*! Ca over Si abundance divided by the solar ratio for these elements */
+ float Ca_over_Si_ratio_in_solar;
- /*! S over Si abundance ratio */
- float sulphur_over_silicon_ratio;
+ /*! S over Si abundance divided by the solar ratio for these elements */
+ float S_over_Si_ratio_in_solar;
/*! Redshift of He reionization */
float He_reion_z_centre;
diff --git a/src/cooling/EAGLE/cooling_tables.c b/src/cooling/EAGLE/cooling_tables.c
index 4f3aed05746ee0f7f6559936ad151d8227bbb90c..e5798a6946fcd52240d4c3985dbc7570e3d8709a 100644
--- a/src/cooling/EAGLE/cooling_tables.c
+++ b/src/cooling/EAGLE/cooling_tables.c
@@ -31,6 +31,7 @@
#include <string.h>
/* Local includes. */
+#include "chemistry_struct.h"
#include "cooling_struct.h"
#include "cooling_tables.h"
#include "error.h"
@@ -39,7 +40,7 @@
/**
* @brief Names of the elements in the order they are stored in the files
*/
-static const char *eagle_tables_element_names[9] = {
+static const char *eagle_tables_element_names[eagle_cooling_N_metal] = {
"Carbon", "Nitrogen", "Oxygen", "Neon", "Magnesium",
"Silicon", "Sulphur", "Calcium", "Iron"};
@@ -211,6 +212,10 @@ void read_cooling_header(const char *fname,
if (N_SolarAbundances != eagle_cooling_N_abundances)
error("Invalid solar abundances array length.");
+ /* Check value */
+ if (N_SolarAbundances != chemistry_element_count + 2)
+ error("Number of abundances not compatible with the chemistry model.");
+
dataset = H5Dopen(tempfile_id, "/Header/Number_of_metals", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&N_Elements);
@@ -237,6 +242,10 @@ void read_cooling_header(const char *fname,
if (posix_memalign((void **)&cooling->SolarAbundances, SWIFT_STRUCT_ALIGNMENT,
N_SolarAbundances * sizeof(float)) != 0)
error("Failed to allocate Solar abundances table");
+ if (posix_memalign((void **)&cooling->SolarAbundances_inv,
+ SWIFT_STRUCT_ALIGNMENT,
+ N_SolarAbundances * sizeof(float)) != 0)
+ error("Failed to allocate Solar abundances inverses table");
/* read in values for each of the arrays */
dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
@@ -285,6 +294,10 @@ void read_cooling_header(const char *fname,
for (int i = 0; i < N_nH; i++) {
cooling->nH[i] = log10(cooling->nH[i]);
}
+ /* Compute inverse of solar mass fractions */
+ for (int i = 0; i < N_SolarAbundances; ++i) {
+ cooling->SolarAbundances_inv[i] = 1.f / cooling->SolarAbundances[i];
+ }
#else
error("Need HDF5 to read cooling tables");
diff --git a/src/scheduler.c b/src/scheduler.c
index 2e56c8b4b0f12c67de54f2fe6043b281c7681523..4e3eb4e29e6cd4a2cd91032d4ee81d203977a59e 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -183,10 +183,11 @@ struct task_dependency {
};
#ifdef WITH_MPI
+
/**
* @brief Define the #task_dependency for MPI
*
- * @param tstype The #MPI_Datatype to initialize
+ * @param tstype The MPI_Datatype to initialize
*/
void task_dependency_define(MPI_Datatype *tstype) {
diff --git a/src/task.c b/src/task.c
index 0ac912b300a3a72e53ac0552ebdec378bb786f52..f16aadc8afb7a2f811c4790688fb849ba1601ce3 100644
--- a/src/task.c
+++ b/src/task.c
@@ -630,7 +630,7 @@ void task_print(const struct task *t) {
* graph.
*
* @param type The #task type.
- * @param subtype The #subtask type.
+ * @param subtype The #task subtype.
* @param cluster (return) The group name (should be allocated)
*/
void task_get_group_name(int type, int subtype, char *cluster) {
diff --git a/src/tracers/EAGLE/tracers_io.h b/src/tracers/EAGLE/tracers_io.h
index 8e662f24ffc4bd0865980975ff5fe1c4f7ddc9a0..492997f53564707a9d8da4dcc3aecefce0edf4ea 100644
--- a/src/tracers/EAGLE/tracers_io.h
+++ b/src/tracers/EAGLE/tracers_io.h
@@ -32,7 +32,6 @@
* @brief Writes the current model of tracers to the file.
*
* @param h_grp The HDF5 group in which to write
- * @param tracers The #tracers_function_data
*/
__attribute__((always_inline)) INLINE static void tracers_write_flavour(
hid_t h_grp) {
@@ -47,6 +46,7 @@ __attribute__((always_inline)) INLINE static void tracers_write_flavour(
* @param parts The particle array.
* @param xparts The extended data particle array.
* @param list The list of i/o properties to write.
+ * @param with_cosmology Are we running with cosmology switched on?
*
* @return Returns the number of fields to write.
*/
diff --git a/tools/analyse_runtime.py b/tools/analyse_runtime.py
index e66fcb0e419b0ae8ae119a4e5f76ad14ae6309a7..f2f198dfb80d6373e63296b6350fe6768191dd39 100755
--- a/tools/analyse_runtime.py
+++ b/tools/analyse_runtime.py
@@ -243,7 +243,7 @@ for i in range(len(labels)):
important_times[0] += times[i]
important_ratios[0] += time_ratios[i]
- print(" - '%-40s' (%5d calls): %.4f%%" % (labels[i], counts[i], time_ratios[i] * 100))
+ print(" - '%-40s' (%5d calls, time: %.4fs): %.4f%%" % (labels[i], counts[i], times[i], time_ratios[i] * 100))
# Anything unaccounted for?
print(