From f53a30d7e0d8b02821a3083be8515008eae60b10 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 12:14:52 +0000 Subject: [PATCH 01/22] Clarified calculation of delta_z in the cosmology theory document. --- theory/Cosmology/operators.tex | 28 +++++++++++++++------------- theory/Cosmology/timesteps.tex | 4 ++++ 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/theory/Cosmology/operators.tex b/theory/Cosmology/operators.tex index 89aa32bae..4ec4ea4b5 100644 --- a/theory/Cosmology/operators.tex +++ b/theory/Cosmology/operators.tex @@ -37,19 +37,21 @@ time using $\Delta t_{\rm kick,A} = \Delta t_{\rm operator. They then use $\int H dt$ as the operator, which integrates out trivially. This slight inconsistency with the rest of the time-integration operators is unlikely to lead to any practical - difference.}, whilst the change in energy due to the expansion of -the Universe (first term in eq.~\ref{eq:cosmo_eom_u}) can be computed -using -\begin{equation} - \int_{a_n}^{a_{n+1}} H dt = \int_{a_n}^{a_{n+1}} \frac{da}{a} = - \log{a_{n+1}} - \log{a_n}. -\end{equation} + difference.}. We additionally compute a few other terms +appearing in some viscosity terms and subgrid models. There are the +difference in cosmic time between the start and the end of the step +and the corresponding change in redshift: +\begin{align} + \Delta t_{\rm cosmic} &= \int_{a_n}^{a_{n+1}} dt = \frac{1}{H_0} + \int_{a_n}^{a_{n+1}} \frac{da}{a E(a)},\\ + \Delta z &= \frac{1}{a_n} - \frac{1}{a_{n+1}} \approx -\frac{H}{a} \Delta t_{\rm cosmic}. +\end{align} Following the same method as for the age of the Universe -(sec. \ref{ssec:flrw}), the three non-trivial integrals are evaluated -numerically at the start of the simulation for a series $10^4$ values -of $a$ placed at regular intervals between $\log a_{\rm begin}$ and -$\log a_{\rm end}$. The values for a specific pair of scale-factors -$a_n$ and $a_{n+1}$ are then obtained by interpolating that table -linearly. +(sec. \ref{ssec:flrw}), these three non-trivial integrals are +evaluated numerically at the start of the simulation for a series +$10^4$ values of $a$ placed at regular intervals between $\log a_{\rm + begin}$ and $\log a_{\rm end}$. The values for a specific pair of +scale-factors $a_n$ and $a_{n+1}$ are then obtained by interpolating +that table linearly. diff --git a/theory/Cosmology/timesteps.tex b/theory/Cosmology/timesteps.tex index 4a1c2ef53..a9856d03a 100644 --- a/theory/Cosmology/timesteps.tex +++ b/theory/Cosmology/timesteps.tex @@ -40,3 +40,7 @@ dominates the overall time-step size calculation. \subsubsection{Conversion from time to integer time-line} +\begin{equation} + \int_{a_n}^{a_{n+1}} H dt = \int_{a_n}^{a_{n+1}} \frac{da}{a} = + \log{a_{n+1}} - \log{a_n}. +\end{equation} -- GitLab From 0525c736a336ea3ea20ee6e39b7e954193ce8053 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 12:17:07 +0000 Subject: [PATCH 02/22] Clarified expression for delta_z in the cooling routine. Added function stubs for temperature calculation. --- src/cooling/EAGLE/cooling.c | 38 ++++++++++++++++++++++++++++--------- src/cooling/EAGLE/cooling.h | 7 +++++++ 2 files changed, 36 insertions(+), 9 deletions(-) diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c index 8dcef4035..1f27d9ca6 100644 --- a/src/cooling/EAGLE/cooling.c +++ b/src/cooling/EAGLE/cooling.c @@ -457,8 +457,11 @@ void cooling_cool_part(const struct phys_const *restrict phys_const, float abundance_ratio[chemistry_element_count + 2]; abundance_ratio_to_solar(p, cooling, abundance_ratio); - /* Get the H and He mass fractions */ + /* Get the Hydrogen mass fraction */ const float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H]; + + /* Get the Helium mass fraction. Note that this is He / (H + He), i.e. a + * metal-free Helium mass fraction as per the Wiersma+08 definition */ const float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] / (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]); @@ -472,14 +475,6 @@ void cooling_cool_part(const struct phys_const *restrict phys_const, * equivalent expression below */ const double ratefact_cgs = n_H_cgs * (XH * cooling->inv_proton_mass_cgs); - /* Get helium and hydrogen reheating term */ - const double Helium_reion_heat_cgs = eagle_helium_reionization_extraheat( - cooling->z_index, -dt * cosmo->H * cosmo->a_inv, cooling); - - /* Convert this into a rate */ - const double Lambda_He_reion_cgs = - Helium_reion_heat_cgs / (dt_cgs * ratefact_cgs); - /* compute hydrogen number density and helium fraction table indices and * offsets (These are fixed for of u, so no need to recompute them) */ int He_i, n_h_i; @@ -488,6 +483,21 @@ void cooling_cool_part(const struct phys_const *restrict phys_const, get_index_1d(cooling->nH, eagle_cooling_N_density, log10(n_H_cgs), &n_h_i, &d_n_h); + /* Start by computing the cooling (heating actually) rate from Helium + re-ionization as this needs to be added on no matter what */ + + /* Change in redshift over the course of this time-step + (See cosmology theory document for the derivation) */ + const double delta_z = -dt * cosmo->H * cosmo->a_inv; + + /* Get helium and hydrogen reheating term */ + const double Helium_reion_heat_cgs = + eagle_helium_reionization_extraheat(cooling->z_index, delta_z, cooling); + + /* Convert this into a rate */ + const double Lambda_He_reion_cgs = + Helium_reion_heat_cgs / (dt_cgs * ratefact_cgs); + /* Let's compute the internal energy at the end of the step */ double u_final_cgs; @@ -620,6 +630,16 @@ __attribute__((always_inline)) INLINE void cooling_first_init_part( xp->cooling_data.radiated_energy = 0.f; } +float cooling_get_temperature( + const struct phys_const *restrict phys_const, + const struct unit_system *restrict us, + const struct cosmology *restrict cosmo, + const struct cooling_function_data *restrict cooling, + const struct part *restrict p, struct xpart *restrict xp) { + + return 1.f; +} + /** * @brief Returns the total radiated energy by this particle. * diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h index 5685692c3..7485880a1 100644 --- a/src/cooling/EAGLE/cooling.h +++ b/src/cooling/EAGLE/cooling.h @@ -59,6 +59,13 @@ void cooling_first_init_part( const struct cooling_function_data *restrict cooling, const struct part *restrict p, struct xpart *restrict xp); +float cooling_get_temperature( + const struct phys_const *restrict phys_const, + const struct unit_system *restrict us, + const struct cosmology *restrict cosmo, + const struct cooling_function_data *restrict cooling, + const struct part *restrict p, struct xpart *restrict xp); + float cooling_get_radiated_energy(const struct xpart *restrict xp); void cooling_init_backend(struct swift_params *parameter_file, -- GitLab From 1411c77055f70e3a323e7083dfdf77feea1bbc54 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 12:35:16 +0000 Subject: [PATCH 03/22] Slight re-organisation of terms inside the main EAGLE cooling function. --- src/cooling/EAGLE/cooling.c | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c index 1f27d9ca6..7325a4bf9 100644 --- a/src/cooling/EAGLE/cooling.c +++ b/src/cooling/EAGLE/cooling.c @@ -453,6 +453,10 @@ void cooling_cool_part(const struct phys_const *restrict phys_const, const double u_0_cgs = u_0 * cooling->internal_energy_to_cgs; const double dt_cgs = dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME); + /* Change in redshift over the course of this time-step + (See cosmology theory document for the derivation) */ + const double delta_redshift = -dt * cosmo->H * cosmo->a_inv; + /* Get this particle's abundance ratios */ float abundance_ratio[chemistry_element_count + 2]; abundance_ratio_to_solar(p, cooling, abundance_ratio); @@ -476,7 +480,8 @@ void cooling_cool_part(const struct phys_const *restrict phys_const, const double ratefact_cgs = n_H_cgs * (XH * cooling->inv_proton_mass_cgs); /* compute hydrogen number density and helium fraction table indices and - * offsets (These are fixed for of u, so no need to recompute them) */ + * offsets (These are fixed for any value of u, so no need to recompute them) + */ int He_i, n_h_i; float d_He, d_n_h; get_index_1d(cooling->HeFrac, eagle_cooling_N_He_frac, HeFrac, &He_i, &d_He); @@ -486,13 +491,9 @@ void cooling_cool_part(const struct phys_const *restrict phys_const, /* Start by computing the cooling (heating actually) rate from Helium re-ionization as this needs to be added on no matter what */ - /* Change in redshift over the course of this time-step - (See cosmology theory document for the derivation) */ - const double delta_z = -dt * cosmo->H * cosmo->a_inv; - /* Get helium and hydrogen reheating term */ const double Helium_reion_heat_cgs = - eagle_helium_reionization_extraheat(cooling->z_index, delta_z, cooling); + eagle_helium_reionization_extraheat(cooling->z_index, delta_redshift, cooling); /* Convert this into a rate */ const double Lambda_He_reion_cgs = -- GitLab From f2b2bffcb09d107b7197c9761556cbafd88df7cb Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 14:45:11 +0000 Subject: [PATCH 04/22] Added functions to compute 10^x even in the absence of GNU extensions. --- src/Makefile.am | 2 +- src/pow10.h | 63 +++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 64 insertions(+), 1 deletion(-) create mode 100644 src/pow10.h diff --git a/src/Makefile.am b/src/Makefile.am index 78531f1f0..7c5ac76bd 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -48,7 +48,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \ dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \ gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \ chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \ - mesh_gravity.h cbrt.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h \ + mesh_gravity.h cbrt.h pow10.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h \ logger_io.h # source files for EAGLE cooling diff --git a/src/pow10.h b/src/pow10.h new file mode 100644 index 000000000..02a9b33b2 --- /dev/null +++ b/src/pow10.h @@ -0,0 +1,63 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@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 . + * + ******************************************************************************/ +#ifndef SWIFT_POW10_H +#define SWIFT_POW10_H + +/* Config parameters. */ +#include "../config.h" + +/* Some standard headers. */ +#include + +/* Local headers. */ +#include "inline.h" + +#ifndef __GNUC__ + +/** + * @brief Raises 10 to the power of the argument. + * + * This function is only used as a replacement for compilers that do + * not implement GNU extensions to the C language. + * + * @param x The input value. + */ +__attribute__((always_inline, const)) INLINE static double pow10( + const double x) { + + return exp(x * M_LN10); +} + +/** + * @brief Raises 10 to the power of the argument. + * + * This function is only used as a replacement for compilers that do + * not implement GNU extensions to the C language. + * + * @param x The input value. + */ +__attribute__((always_inline, const)) INLINE static float pow10f( + const float x) { + + return expf(x * (float)M_LN10); +} + +#endif /* __GNUC__ */ + +#endif /* SWIFT_POW10_H */ -- GitLab From 77e3b8064e3f2b54e753da66047953848c5f835d Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 14:53:37 +0000 Subject: [PATCH 05/22] The function should be called exp10 and now pow10 --- src/Makefile.am | 2 +- src/exp10.h | 63 +++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 64 insertions(+), 1 deletion(-) create mode 100644 src/exp10.h diff --git a/src/Makefile.am b/src/Makefile.am index 7c5ac76bd..889b0a3af 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -48,7 +48,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \ dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \ gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \ chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \ - mesh_gravity.h cbrt.h pow10.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h \ + mesh_gravity.h cbrt.h exp10.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h \ logger_io.h # source files for EAGLE cooling diff --git a/src/exp10.h b/src/exp10.h new file mode 100644 index 000000000..48a950574 --- /dev/null +++ b/src/exp10.h @@ -0,0 +1,63 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@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 . + * + ******************************************************************************/ +#ifndef SWIFT_EXP10_H +#define SWIFT_EXP10_H + +/* Config parameters. */ +#include "../config.h" + +/* Some standard headers. */ +#include + +/* Local headers. */ +#include "inline.h" + +#ifndef __GNUC__ + +/** + * @brief Raises 10 to the power of the argument. + * + * This function is only used as a replacement for compilers that do + * not implement GNU extensions to the C language. + * + * @param x The input value. + */ +__attribute__((always_inline, const)) INLINE static double exp10( + const double x) { + + return exp(x * M_LN10); +} + +/** + * @brief Raises 10 to the power of the argument. + * + * This function is only used as a replacement for compilers that do + * not implement GNU extensions to the C language. + * + * @param x The input value. + */ +__attribute__((always_inline, const)) INLINE static float exp10f( + const float x) { + + return expf(x * (float)M_LN10); +} + +#endif /* __GNUC__ */ + +#endif /* SWIFT_EXP10_H */ -- GitLab From f3abd9a56fc9ed03d177a135442d88c7a7f4bf07 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 14:54:13 +0000 Subject: [PATCH 06/22] The function should be called exp10 and now pow10 --- src/pow10.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pow10.h b/src/pow10.h index 02a9b33b2..3414072a8 100644 --- a/src/pow10.h +++ b/src/pow10.h @@ -38,7 +38,7 @@ * * @param x The input value. */ -__attribute__((always_inline, const)) INLINE static double pow10( +__attribute__((always_inline, const)) INLINE static double exp10( const double x) { return exp(x * M_LN10); @@ -52,7 +52,7 @@ __attribute__((always_inline, const)) INLINE static double pow10( * * @param x The input value. */ -__attribute__((always_inline, const)) INLINE static float pow10f( +__attribute__((always_inline, const)) INLINE static float exp10f( const float x) { return expf(x * (float)M_LN10); -- GitLab From 8e2729b1352c8f72b13776589e0ef133f1e31d66 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 15:00:06 +0000 Subject: [PATCH 07/22] Use the Temperature and not log10(T) when computing Compton cooling. --- src/cooling/EAGLE/cooling_rates.h | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/cooling/EAGLE/cooling_rates.h b/src/cooling/EAGLE/cooling_rates.h index 621f7c7b0..c9de0d642 100644 --- a/src/cooling/EAGLE/cooling_rates.h +++ b/src/cooling/EAGLE/cooling_rates.h @@ -25,6 +25,7 @@ /* Local includes. */ #include "cooling_tables.h" +#include "exp10.h" #include "interpolate.h" /** @@ -364,14 +365,15 @@ INLINE static double eagle_metal_cooling_rate( /* Temperature */ float dT_du = -1.f; - const double T = + const double log_10_T = eagle_convert_u_to_temp(log10_u_cgs, redshift, compute_dT_du, &dT_du, n_H_index, He_index, d_n_h, d_He, cooling); /* Get index along temperature dimension of the tables */ int T_index; float d_T; - get_index_1d(cooling->Temp, eagle_cooling_N_temperature, T, &T_index, &d_T); + get_index_1d(cooling->Temp, eagle_cooling_N_temperature, log_10_T, &T_index, + &d_T); #ifdef TO_BE_DONE /* Difference between entries on the temperature table around u */ @@ -516,6 +518,8 @@ INLINE static double eagle_metal_cooling_rate( if ((redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) || (redshift > cooling->H_reion_z)) { + const double T = exp10(log_10_T); + /* Note the minus sign */ Lambda_Compton -= eagle_Compton_cooling_rate(cooling, redshift, n_H_cgs, T, H_plus_He_electron_abundance); -- GitLab From 9838e226ddec0b774dbdb682c6c35b513f7a3f9d Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 15:06:50 +0000 Subject: [PATCH 08/22] Added function to calculate temperature from the Wiersma table interpolation. --- src/cooling/EAGLE/cooling.c | 39 ++++++++++++++++++++++++++++++++++--- 1 file changed, 36 insertions(+), 3 deletions(-) diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c index 7325a4bf9..d5bb60d58 100644 --- a/src/cooling/EAGLE/cooling.c +++ b/src/cooling/EAGLE/cooling.c @@ -492,8 +492,8 @@ void cooling_cool_part(const struct phys_const *restrict phys_const, re-ionization as this needs to be added on no matter what */ /* Get helium and hydrogen reheating term */ - const double Helium_reion_heat_cgs = - eagle_helium_reionization_extraheat(cooling->z_index, delta_redshift, cooling); + const double Helium_reion_heat_cgs = eagle_helium_reionization_extraheat( + cooling->z_index, delta_redshift, cooling); /* Convert this into a rate */ const double Lambda_He_reion_cgs = @@ -638,7 +638,40 @@ float cooling_get_temperature( const struct cooling_function_data *restrict cooling, const struct part *restrict p, struct xpart *restrict xp) { - return 1.f; + /* Get physical internal energy */ + const float u = hydro_get_physical_internal_energy(p, xp, cosmo); + const double u_cgs = u * cooling->internal_energy_to_cgs; + + /* Get the Hydrogen mass fraction */ + const float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H]; + + /* Get the Helium mass fraction. Note that this is He / (H + He), i.e. a + * metal-free Helium mass fraction as per the Wiersma+08 definition */ + const float HeFrac = + p->chemistry_data.metal_mass_fraction[chemistry_element_He] / + (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]); + + /* Convert Hydrogen mass fraction into Hydrogen number density */ + const float rho = hydro_get_physical_density(p, cosmo); + const double n_H = rho * XH / phys_const->const_proton_mass; + const double n_H_cgs = n_H * cooling->number_density_to_cgs; + + /* compute hydrogen number density and helium fraction table indices and + * offsets */ + int He_index, n_H_index; + float d_He, d_n_H; + get_index_1d(cooling->HeFrac, eagle_cooling_N_He_frac, HeFrac, &He_index, + &d_He); + get_index_1d(cooling->nH, eagle_cooling_N_density, log10(n_H_cgs), &n_H_index, + &d_n_H); + + /* Compute the log10 of the temperature by interpolating the table */ + const double log_10_T = eagle_convert_u_to_temp( + log10(u_cgs), cosmo->z, /*compute_dT_du=*/0, /*dT_du=*/NULL, n_H_index, + He_index, d_n_H, d_He, cooling); + + /* Undo the log! */ + return exp10(log_10_T); } /** -- GitLab From 46ea4e845fc383d1d017f1aefcf32f4018d1aaef Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 15:10:26 +0000 Subject: [PATCH 09/22] Removed unnecessarily added file pow10.h --- src/pow10.h | 63 ----------------------------------------------------- 1 file changed, 63 deletions(-) delete mode 100644 src/pow10.h diff --git a/src/pow10.h b/src/pow10.h deleted file mode 100644 index 3414072a8..000000000 --- a/src/pow10.h +++ /dev/null @@ -1,63 +0,0 @@ -/******************************************************************************* - * This file is part of SWIFT. - * Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@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 . - * - ******************************************************************************/ -#ifndef SWIFT_POW10_H -#define SWIFT_POW10_H - -/* Config parameters. */ -#include "../config.h" - -/* Some standard headers. */ -#include - -/* Local headers. */ -#include "inline.h" - -#ifndef __GNUC__ - -/** - * @brief Raises 10 to the power of the argument. - * - * This function is only used as a replacement for compilers that do - * not implement GNU extensions to the C language. - * - * @param x The input value. - */ -__attribute__((always_inline, const)) INLINE static double exp10( - const double x) { - - return exp(x * M_LN10); -} - -/** - * @brief Raises 10 to the power of the argument. - * - * This function is only used as a replacement for compilers that do - * not implement GNU extensions to the C language. - * - * @param x The input value. - */ -__attribute__((always_inline, const)) INLINE static float exp10f( - const float x) { - - return expf(x * (float)M_LN10); -} - -#endif /* __GNUC__ */ - -#endif /* SWIFT_POW10_H */ -- GitLab From f2f970c45609defbf98fd78f1961bcccb96d7c0d Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 15:18:35 +0000 Subject: [PATCH 10/22] Added documentation of the function returning temperatures in the EAGLE model. --- src/cooling/EAGLE/cooling.c | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c index d5bb60d58..5324ae47d 100644 --- a/src/cooling/EAGLE/cooling.c +++ b/src/cooling/EAGLE/cooling.c @@ -37,6 +37,7 @@ #include "cooling_struct.h" #include "cooling_tables.h" #include "error.h" +#include "exp10.h" #include "hydro.h" #include "interpolate.h" #include "io_properties.h" @@ -631,6 +632,22 @@ __attribute__((always_inline)) INLINE void cooling_first_init_part( xp->cooling_data.radiated_energy = 0.f; } +/** + * @brief Compute the temperature of a #part based on the cooling function. + * + * We use the Temperature table of the Wiersma+08 set. This computes the + * equilibirum temperature of a gas for a given redshift, Hydrogen density, + * internal energy per unit mass and Helium fraction. + * + * The temperature returned is consistent with the cooling rates. + * + * @param phys_const #phys_const data structure. + * @param us The internal system of units. + * @param cosmo #cosmology data structure. + * @param cooling #cooling_function_data struct. + * @param p #part data. + * @param xp Pointer to the #xpart data. + */ float cooling_get_temperature( const struct phys_const *restrict phys_const, const struct unit_system *restrict us, -- GitLab From 01a56b9dd344722dbdf2f10facb689c1665d5be2 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 15:35:20 +0000 Subject: [PATCH 11/22] Updated all the YAML files with the correct parameter names for the EAGLE cooling model. --- examples/EAGLE_12/eagle_12.yml | 9 ++++----- examples/EAGLE_25/eagle_25.yml | 8 ++++---- examples/EAGLE_50/eagle_50.yml | 8 ++++---- examples/EAGLE_6/eagle_6.yml | 8 ++++---- examples/parameter_example.yml | 7 +++++++ 5 files changed, 23 insertions(+), 17 deletions(-) diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml index 90b546f31..1d1932049 100644 --- a/examples/EAGLE_12/eagle_12.yml +++ b/examples/EAGLE_12/eagle_12.yml @@ -72,9 +72,8 @@ EAGLEChemistry: # Solar abundances SulphurOverSilicon: 0.6054160 EagleCooling: - filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/ + filename: ./coolingtables/ reionisation_redshift: 11.5 - he_reion_z_center: 3.5 - he_reion_z_sigma: 0.5 - he_reion_ev_pH: 2.0 - + He_reion_z_centre: 3.5 + He_reion_z_sigma: 0.5 + He_reion_ev_pH: 2.0 diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml index bd74473d1..3e1e5907d 100644 --- a/examples/EAGLE_25/eagle_25.yml +++ b/examples/EAGLE_25/eagle_25.yml @@ -81,9 +81,9 @@ EAGLEChemistry: # Solar abundances SulphurOverSilicon: 0.6054160 EagleCooling: - filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/ + filename: ./coolingtables/ reionisation_redshift: 11.5 - he_reion_z_center: 3.5 - he_reion_z_sigma: 0.5 - he_reion_ev_pH: 2.0 + He_reion_z_centre: 3.5 + He_reion_z_sigma: 0.5 + He_reion_ev_pH: 2.0 diff --git a/examples/EAGLE_50/eagle_50.yml b/examples/EAGLE_50/eagle_50.yml index b86a3d87d..87bef197c 100644 --- a/examples/EAGLE_50/eagle_50.yml +++ b/examples/EAGLE_50/eagle_50.yml @@ -74,8 +74,8 @@ EAGLEChemistry: # Solar abundances SulphurOverSilicon: 0.6054160 EagleCooling: - filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/ + filename: ./coolingtables/ reionisation_redshift: 11.5 - he_reion_z_center: 3.5 - he_reion_z_sigma: 0.5 - he_reion_ev_pH: 2.0 + He_reion_z_centre: 3.5 + He_reion_z_sigma: 0.5 + He_reion_ev_pH: 2.0 \ No newline at end of file diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml index 494f48b83..95a5d3398 100644 --- a/examples/EAGLE_6/eagle_6.yml +++ b/examples/EAGLE_6/eagle_6.yml @@ -85,8 +85,8 @@ EAGLEChemistry: # Solar abundances SulphurOverSilicon: 0.6054160 EagleCooling: - filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/ + filename: ./coolingtables/ reionisation_redshift: 11.5 - he_reion_z_center: 3.5 - he_reion_z_sigma: 0.5 - he_reion_ev_pH: 2.0 + He_reion_z_centre: 3.5 + He_reion_z_sigma: 0.5 + He_reion_ev_pH: 2.0 \ No newline at end of file diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 6adccf296..d7b708d54 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -245,6 +245,13 @@ 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. +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 + # Cooling with Grackle 3.0 GrackleCooling: CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository) -- GitLab From 645ffcdfd99edccdc6a9ca6cc3bbe8449bb145e6 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 15:43:34 +0000 Subject: [PATCH 12/22] The function computing the gas temperature should not modify the xpart. --- src/cooling/EAGLE/cooling.c | 2 +- src/cooling/EAGLE/cooling.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c index 5324ae47d..0815b6cb8 100644 --- a/src/cooling/EAGLE/cooling.c +++ b/src/cooling/EAGLE/cooling.c @@ -653,7 +653,7 @@ float cooling_get_temperature( const struct unit_system *restrict us, const struct cosmology *restrict cosmo, const struct cooling_function_data *restrict cooling, - const struct part *restrict p, struct xpart *restrict xp) { + const struct part *restrict p, const struct xpart *restrict xp) { /* Get physical internal energy */ const float u = hydro_get_physical_internal_energy(p, xp, cosmo); diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h index 7485880a1..733749548 100644 --- a/src/cooling/EAGLE/cooling.h +++ b/src/cooling/EAGLE/cooling.h @@ -64,7 +64,7 @@ float cooling_get_temperature( const struct unit_system *restrict us, const struct cosmology *restrict cosmo, const struct cooling_function_data *restrict cooling, - const struct part *restrict p, struct xpart *restrict xp); + const struct part *restrict p, const struct xpart *restrict xp); float cooling_get_radiated_energy(const struct xpart *restrict xp); -- GitLab From 9b4f1c0f1c431ec952ea03d7072201ca464dfdc7 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 15:45:13 +0000 Subject: [PATCH 13/22] Update the small cosmo volume YAML files with the correct parameter names for the EAGLE cooling model. --- examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml b/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml index 8ad9ae074..969626d97 100644 --- a/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml +++ b/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml @@ -66,10 +66,9 @@ LambdaCooling: EagleCooling: filename: ./coolingtables/ reionisation_redshift: 8.5 - he_reion: 1 - he_reion_z_center: 3.5 - he_reion_z_sigma: 0.5 - he_reion_ev_pH: 2.0 + He_reion_z_centre: 3.5 + He_reion_z_sigma: 0.5 + He_reion_ev_pH: 2.0 # Impose primoridal metallicity EAGLEChemistry: -- GitLab From 8586ddbb866a00bb32d9d2852bb3a62c889bd983 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Mon, 10 Dec 2018 16:10:31 +0000 Subject: [PATCH 14/22] Add the Temperature to the output of the particles in the EAGLE cooling model. --- src/cooling/EAGLE/cooling_io.h | 17 +++++++++++++++-- src/parallel_io.c | 9 +++++---- src/serial_io.c | 9 +++++---- src/single_io.c | 9 +++++---- 4 files changed, 30 insertions(+), 14 deletions(-) diff --git a/src/cooling/EAGLE/cooling_io.h b/src/cooling/EAGLE/cooling_io.h index 48c845c25..86b7067bd 100644 --- a/src/cooling/EAGLE/cooling_io.h +++ b/src/cooling/EAGLE/cooling_io.h @@ -23,6 +23,7 @@ #include "../config.h" /* Local includes */ +#include "cooling.h" #include "io_properties.h" #ifdef HAVE_HDF5 @@ -40,9 +41,17 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( } #endif +INLINE static void convert_part_T(const struct engine* e, const struct part* p, + const struct xpart* xp, float* ret) { + + ret[0] = cooling_get_temperature(e->physical_constants, e->internal_units, + e->cosmology, e->cooling_func, p, xp); +} + /** * @brief Specifies which particle fields to write to a dataset * + * @param parts The particle array. * @param xparts The extended data particle array. * @param list The list of i/o properties to write. * @param cooling The #cooling_function_data @@ -50,9 +59,13 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( * @return Returns the number of fields to write. */ __attribute__((always_inline)) INLINE static int cooling_write_particles( - const struct xpart* xparts, struct io_props* list, + const struct part* parts, const struct xpart* xparts, struct io_props* list, const struct cooling_function_data* cooling) { - return 0; + + list[0] = io_make_output_field_convert_part("Temperature", FLOAT, 1, + UNIT_CONV_TEMPERATURE, parts, + xparts, convert_part_T); + return 1; } #endif /* SWIFT_COOLING_EAGLE_IO_H */ diff --git a/src/parallel_io.c b/src/parallel_io.c index febbb4a7d..5e616d52b 100644 --- a/src/parallel_io.c +++ b/src/parallel_io.c @@ -1394,8 +1394,8 @@ void write_output_parallel(struct engine* e, const char* baseName, Nparticles = Ngas; hydro_write_particles(parts, xparts, list, &num_fields); num_fields += chemistry_write_particles(parts, list + num_fields); - num_fields += cooling_write_particles(xparts, list + num_fields, - e->cooling_func); + num_fields += cooling_write_particles( + parts, xparts, list + num_fields, e->cooling_func); } else { /* Ok, we need to fish out the particles we want */ @@ -1418,8 +1418,9 @@ void write_output_parallel(struct engine* e, const char* baseName, &num_fields); num_fields += chemistry_write_particles(parts_written, list + num_fields); - num_fields += cooling_write_particles( - xparts_written, list + num_fields, e->cooling_func); + num_fields += + cooling_write_particles(parts_written, xparts_written, + list + num_fields, e->cooling_func); } } break; diff --git a/src/serial_io.c b/src/serial_io.c index 059318df1..76f55897d 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -1071,8 +1071,8 @@ void write_output_serial(struct engine* e, const char* baseName, Nparticles = Ngas; hydro_write_particles(parts, xparts, list, &num_fields); num_fields += chemistry_write_particles(parts, list + num_fields); - num_fields += cooling_write_particles(xparts, list + num_fields, - e->cooling_func); + num_fields += cooling_write_particles( + parts, xparts, list + num_fields, e->cooling_func); } else { /* Ok, we need to fish out the particles we want */ @@ -1095,8 +1095,9 @@ void write_output_serial(struct engine* e, const char* baseName, &num_fields); num_fields += chemistry_write_particles(parts_written, list + num_fields); - num_fields += cooling_write_particles( - xparts_written, list + num_fields, e->cooling_func); + num_fields += + cooling_write_particles(parts_written, xparts_written, + list + num_fields, e->cooling_func); } } break; diff --git a/src/single_io.c b/src/single_io.c index 833c4a80c..67c3250e3 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -890,8 +890,8 @@ void write_output_single(struct engine* e, const char* baseName, N = Ngas; hydro_write_particles(parts, xparts, list, &num_fields); num_fields += chemistry_write_particles(parts, list + num_fields); - num_fields += cooling_write_particles(xparts, list + num_fields, - e->cooling_func); + num_fields += cooling_write_particles( + parts, xparts, list + num_fields, e->cooling_func); } else { /* Ok, we need to fish out the particles we want */ @@ -914,8 +914,9 @@ void write_output_single(struct engine* e, const char* baseName, &num_fields); num_fields += chemistry_write_particles(parts_written, list + num_fields); - num_fields += cooling_write_particles( - xparts_written, list + num_fields, e->cooling_func); + num_fields += + cooling_write_particles(parts_written, xparts_written, + list + num_fields, e->cooling_func); } } break; -- GitLab From 4686c6ec254c017f492305a98bf7a7d635d40879 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Tue, 11 Dec 2018 09:20:52 +0000 Subject: [PATCH 15/22] d_n_h --> d_n_H and n_h_i --> n_H_index in the EAGLE cooling model. --- src/cooling/EAGLE/cooling.c | 82 ++++++++++++++++--------------- src/cooling/EAGLE/cooling_rates.h | 28 +++++------ 2 files changed, 56 insertions(+), 54 deletions(-) diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c index 0815b6cb8..cd3d80efb 100644 --- a/src/cooling/EAGLE/cooling.c +++ b/src/cooling/EAGLE/cooling.c @@ -152,9 +152,9 @@ void cooling_update(const struct cosmology *cosmo, * * @param logu_init Initial guess for log(internal energy) * @param u_ini Internal energy at beginning of hydro step - * @param n_h_i Particle hydrogen number density index - * @param d_n_h Particle hydrogen number density offset - * @param He_i Particle helium fraction index + * @param n_H_index Particle hydrogen number density index + * @param d_n_H Particle hydrogen number density offset + * @param He_index Particle helium fraction index * @param d_He Particle helium fraction offset * @param He_reion_heat Heating due to helium reionization * (only depends on redshift, so passed as parameter) @@ -167,8 +167,8 @@ void cooling_update(const struct cosmology *cosmo, * @param bisection_flag Flag to identify if scheme failed to converge */ INLINE static float newton_iter( - float logu_init, double u_ini, int n_h_i, float d_n_h, int He_i, float d_He, - float He_reion_heat, struct part *restrict p, + float logu_init, double u_ini, int n_H_index, float d_n_H, int He_index, + float d_He, float He_reion_heat, struct part *restrict p, const struct cosmology *restrict cosmo, const struct cooling_function_data *restrict cooling, const struct phys_const *restrict phys_const, @@ -203,10 +203,10 @@ INLINE static float newton_iter( { logu_old = logu; LambdaNet_old = LambdaNet; - LambdaNet = - (He_reion_heat / (dt * ratefact_cgs)) + - eagle_cooling_rate(logu_old, cosmo->z, n_H_cgs, abundance_ratio, n_h_i, - d_n_h, He_i, d_He, cooling, &dLambdaNet_du); + LambdaNet = (He_reion_heat / (dt * ratefact_cgs)) + + eagle_cooling_rate(logu_old, cosmo->z, n_H_cgs, abundance_ratio, + n_H_index, d_n_H, He_index, d_He, cooling, + &dLambdaNet_du); /* Newton iteration. For details on how the cooling equation is integrated * see documentation in theory/Cooling/ */ @@ -244,9 +244,9 @@ INLINE static float newton_iter( * @param u_ini_cgs Internal energy at beginning of hydro step in CGS. * @param n_H_cgs Hydrogen number density in CGS. * @param redshift Current redshift. - * @param n_h_i Particle hydrogen number density index. - * @param d_n_h Particle hydrogen number density offset. - * @param He_i Particle helium fraction index. + * @param n_H_index Particle hydrogen number density index. + * @param d_n_H Particle hydrogen number density offset. + * @param He_index Particle helium fraction index. * @param d_He Particle helium fraction offset. * @param Lambda_He_reion_cgs Cooling rate coming from He reionization. * @param ratefact_cgs Multiplication factor to get a cooling rate. @@ -257,8 +257,9 @@ INLINE static float newton_iter( */ INLINE static double bisection_iter( const double u_ini_cgs, const double n_H_cgs, const double redshift, - int n_h_i, float d_n_h, int He_i, float d_He, double Lambda_He_reion_cgs, - double ratefact_cgs, const struct cooling_function_data *restrict cooling, + int n_H_index, float d_n_H, int He_index, float d_He, + double Lambda_He_reion_cgs, double ratefact_cgs, + const struct cooling_function_data *restrict cooling, const float abundance_ratio[chemistry_element_count + 2], double dt_cgs, long long ID) { @@ -271,10 +272,10 @@ INLINE static double bisection_iter( /*************************************/ double LambdaNet_cgs = - Lambda_He_reion_cgs + eagle_cooling_rate(log(u_ini_cgs), redshift, - n_H_cgs, abundance_ratio, n_h_i, - d_n_h, He_i, d_He, cooling, - /*dLambdaNet_du=*/NULL); + Lambda_He_reion_cgs + + eagle_cooling_rate(log(u_ini_cgs), redshift, n_H_cgs, abundance_ratio, + n_H_index, d_n_H, He_index, d_He, cooling, + /*dLambdaNet_du=*/NULL); /*************************************/ /* Let's try to bracket the solution */ @@ -290,7 +291,7 @@ INLINE static double bisection_iter( LambdaNet_cgs = Lambda_He_reion_cgs + eagle_cooling_rate(log(u_lower_cgs), redshift, n_H_cgs, abundance_ratio, - n_h_i, d_n_h, He_i, d_He, cooling, + n_H_index, d_n_H, He_index, d_He, cooling, /*dLambdaNet_du=*/NULL); int i = 0; @@ -302,11 +303,11 @@ INLINE static double bisection_iter( u_upper_cgs /= bracket_factor; /* Compute a new rate */ - LambdaNet_cgs = - Lambda_He_reion_cgs + - eagle_cooling_rate(log(u_lower_cgs), redshift, n_H_cgs, - abundance_ratio, n_h_i, d_n_h, He_i, d_He, cooling, - /*dLambdaNet_du=*/NULL); + LambdaNet_cgs = Lambda_He_reion_cgs + + eagle_cooling_rate(log(u_lower_cgs), redshift, n_H_cgs, + abundance_ratio, n_H_index, d_n_H, + He_index, d_He, cooling, + /*dLambdaNet_du=*/NULL); i++; } @@ -326,7 +327,7 @@ INLINE static double bisection_iter( LambdaNet_cgs = Lambda_He_reion_cgs + eagle_cooling_rate(log(u_upper_cgs), redshift, n_H_cgs, abundance_ratio, - n_h_i, d_n_h, He_i, d_He, cooling, + n_H_index, d_n_H, He_index, d_He, cooling, /*dLambdaNet_du=*/NULL); int i = 0; @@ -338,11 +339,11 @@ INLINE static double bisection_iter( u_upper_cgs *= bracket_factor; /* Compute a new rate */ - LambdaNet_cgs = - Lambda_He_reion_cgs + - eagle_cooling_rate(log(u_upper_cgs), redshift, n_H_cgs, - abundance_ratio, n_h_i, d_n_h, He_i, d_He, cooling, - /*dLambdaNet_du=*/NULL); + LambdaNet_cgs = Lambda_He_reion_cgs + + eagle_cooling_rate(log(u_upper_cgs), redshift, n_H_cgs, + abundance_ratio, n_H_index, d_n_H, + He_index, d_He, cooling, + /*dLambdaNet_du=*/NULL); i++; } @@ -372,7 +373,7 @@ INLINE static double bisection_iter( LambdaNet_cgs = Lambda_He_reion_cgs + eagle_cooling_rate(log(u_next_cgs), redshift, n_H_cgs, abundance_ratio, - n_h_i, d_n_h, He_i, d_He, cooling, + n_H_index, d_n_H, He_index, d_He, cooling, /*dLambdaNet_du=*/NULL); /* Where do we go next? */ @@ -483,11 +484,12 @@ void cooling_cool_part(const struct phys_const *restrict phys_const, /* compute hydrogen number density and helium fraction table indices and * offsets (These are fixed for any value of u, so no need to recompute them) */ - int He_i, n_h_i; - float d_He, d_n_h; - get_index_1d(cooling->HeFrac, eagle_cooling_N_He_frac, HeFrac, &He_i, &d_He); - get_index_1d(cooling->nH, eagle_cooling_N_density, log10(n_H_cgs), &n_h_i, - &d_n_h); + int He_index, n_H_index; + float d_He, d_n_H; + get_index_1d(cooling->HeFrac, eagle_cooling_N_He_frac, HeFrac, &He_index, + &d_He); + get_index_1d(cooling->nH, eagle_cooling_N_density, log10(n_H_cgs), &n_H_index, + &d_n_H); /* Start by computing the cooling (heating actually) rate from Helium re-ionization as this needs to be added on no matter what */ @@ -506,8 +508,8 @@ void cooling_cool_part(const struct phys_const *restrict phys_const, /* First try an explicit integration (note we ignore the derivative) */ const double LambdaNet_cgs = Lambda_He_reion_cgs + eagle_cooling_rate(log(u_0_cgs), cosmo->z, n_H_cgs, - abundance_ratio, n_h_i, d_n_h, - He_i, d_He, cooling, + abundance_ratio, n_H_index, + d_n_H, He_index, d_He, cooling, /*dLambdaNet_du=*/NULL); /* if cooling rate is small, take the explicit solution */ @@ -545,8 +547,8 @@ void cooling_cool_part(const struct phys_const *restrict phys_const, /* Alright, all else failed, let's bisect */ if (bisection_flag || !(cooling->newton_flag)) { u_final_cgs = - bisection_iter(u_0_cgs, n_H_cgs, cosmo->z, n_h_i, d_n_h, He_i, d_He, - Lambda_He_reion_cgs, ratefact_cgs, cooling, + bisection_iter(u_0_cgs, n_H_cgs, cosmo->z, n_H_index, d_n_H, He_index, + d_He, Lambda_He_reion_cgs, ratefact_cgs, cooling, abundance_ratio, dt_cgs, p->id); } } diff --git a/src/cooling/EAGLE/cooling_rates.h b/src/cooling/EAGLE/cooling_rates.h index c9de0d642..c78be6a30 100644 --- a/src/cooling/EAGLE/cooling_rates.h +++ b/src/cooling/EAGLE/cooling_rates.h @@ -333,7 +333,7 @@ __attribute__((always_inline)) INLINE double eagle_Compton_cooling_rate( * to solar metal abundances * * @param n_H_index Particle hydrogen number density index - * @param d_n_h Particle hydrogen number density offset + * @param d_n_H Particle hydrogen number density offset * @param He_index Particle helium fraction index * @param d_He Particle helium fraction offset * @param cooling Cooling data structure @@ -346,7 +346,7 @@ __attribute__((always_inline)) INLINE double eagle_Compton_cooling_rate( INLINE static double eagle_metal_cooling_rate( double log10_u_cgs, double redshift, double n_H_cgs, const float solar_ratio[chemistry_element_count + 2], int n_H_index, - float d_n_h, int He_index, float d_He, + float d_n_H, int He_index, float d_He, const struct cooling_function_data *restrict cooling, double *dlambda_du, double *element_lambda) { @@ -367,7 +367,7 @@ INLINE static double eagle_metal_cooling_rate( float dT_du = -1.f; const double log_10_T = eagle_convert_u_to_temp(log10_u_cgs, redshift, compute_dT_du, &dT_du, - n_H_index, He_index, d_n_h, d_He, cooling); + n_H_index, He_index, d_n_H, d_He, cooling); /* Get index along temperature dimension of the tables */ int T_index; @@ -393,7 +393,7 @@ INLINE static double eagle_metal_cooling_rate( * in redshift */ Lambda_free = interpolation_3d(cooling->table.H_plus_He_heating, /* */ n_H_index, He_index, T_index, /* */ - d_n_h, d_He, d_T, /* */ + d_n_H, d_He, d_T, /* */ eagle_cooling_N_density, /* */ eagle_cooling_N_He_frac, /* */ eagle_cooling_N_temperature); /* */ @@ -418,7 +418,7 @@ INLINE static double eagle_metal_cooling_rate( Lambda_free = interpolation_4d(cooling->table.H_plus_He_heating, /* */ /*z_index=*/0, n_H_index, He_index, T_index, /* */ - cooling->dz, d_n_h, d_He, d_T, /* */ + cooling->dz, d_n_H, d_He, d_T, /* */ eagle_cooling_N_loaded_redshifts, /* */ eagle_cooling_N_density, /* */ eagle_cooling_N_He_frac, /* */ @@ -462,7 +462,7 @@ INLINE static double eagle_metal_cooling_rate( H_plus_He_electron_abundance = interpolation_3d(cooling->table.H_plus_He_electron_abundance, /* */ n_H_index, He_index, T_index, /* */ - d_n_h, d_He, d_T, /* */ + d_n_H, d_He, d_T, /* */ eagle_cooling_N_density, /* */ eagle_cooling_N_He_frac, /* */ eagle_cooling_N_temperature); /* */ @@ -487,7 +487,7 @@ INLINE static double eagle_metal_cooling_rate( H_plus_He_electron_abundance = interpolation_4d(cooling->table.H_plus_He_electron_abundance, /* */ /*z_index=*/0, n_H_index, He_index, T_index, /* */ - cooling->dz, d_n_h, d_He, d_T, /* */ + cooling->dz, d_n_H, d_He, d_T, /* */ eagle_cooling_N_loaded_redshifts, /* */ eagle_cooling_N_density, /* */ eagle_cooling_N_He_frac, /* */ @@ -543,7 +543,7 @@ INLINE static double eagle_metal_cooling_rate( solar_electron_abundance = interpolation_2d(cooling->table.electron_abundance, /* */ n_H_index, T_index, /* */ - d_n_h, d_T, /* */ + d_n_H, d_T, /* */ eagle_cooling_N_density, /* */ eagle_cooling_N_temperature); /* */ @@ -566,7 +566,7 @@ INLINE static double eagle_metal_cooling_rate( solar_electron_abundance = interpolation_3d(cooling->table.electron_abundance, /* */ /*z_index=*/0, n_H_index, T_index, /* */ - cooling->dz, d_n_h, d_T, /* */ + cooling->dz, d_n_H, d_T, /* */ eagle_cooling_N_loaded_redshifts, /* */ eagle_cooling_N_density, /* */ eagle_cooling_N_temperature); /* */ @@ -605,7 +605,7 @@ INLINE static double eagle_metal_cooling_rate( 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, /* */ + /*delta_elem=*/0.f, d_n_H, d_T, /* */ eagle_cooling_N_metal, /* */ eagle_cooling_N_density, /* */ eagle_cooling_N_temperature); /* */ @@ -641,7 +641,7 @@ INLINE static double eagle_metal_cooling_rate( 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, /* */ + /*delta_elem=*/0.f, cooling->dz, d_n_H, d_T, /* */ eagle_cooling_N_metal, /* */ eagle_cooling_N_loaded_redshifts, /* */ eagle_cooling_N_density, /* */ @@ -700,7 +700,7 @@ INLINE static double eagle_metal_cooling_rate( * @param abundance_ratio Ratio of element abundance to solar. * * @param n_H_index Particle hydrogen number density index - * @param d_n_h Particle hydrogen number density offset + * @param d_n_H Particle hydrogen number density offset * @param He_index Particle helium fraction index * @param d_He Particle helium fraction offset * @param cooling #cooling_function_data structure @@ -713,12 +713,12 @@ INLINE static double eagle_metal_cooling_rate( INLINE static double eagle_cooling_rate( double log_u_cgs, double redshift, double n_H_cgs, const float abundance_ratio[chemistry_element_count + 2], int n_H_index, - float d_n_h, int He_index, float d_He, + float d_n_H, int He_index, float d_He, const struct cooling_function_data *restrict cooling, double *dLambdaNet_du) { return eagle_metal_cooling_rate(log_u_cgs / M_LN10, redshift, n_H_cgs, - abundance_ratio, n_H_index, d_n_h, He_index, + abundance_ratio, n_H_index, d_n_H, He_index, d_He, cooling, dLambdaNet_du, /*element_lambda=*/NULL); } -- GitLab From ecb93a1b77f4c3e2189f83d46ea69d3531d2b161 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Tue, 11 Dec 2018 10:07:04 +0000 Subject: [PATCH 16/22] Add the mean molecular weight of neutral and ionised gas to the hydro properties. --- src/hydro_properties.c | 6 ++++++ src/hydro_properties.h | 6 ++++++ 2 files changed, 12 insertions(+) diff --git a/src/hydro_properties.c b/src/hydro_properties.c index 2b1cd4205..85f88d418 100644 --- a/src/hydro_properties.c +++ b/src/hydro_properties.c @@ -119,6 +119,12 @@ void hydro_props_init(struct hydro_props *p, p->hydrogen_mass_fraction = parser_get_opt_param_double( params, "SPH:H_mass_fraction", default_H_fraction); + /* Mean molecular mass for neutral gas */ + p->mu_neutral = 4. / (1. + 3. * p->hydrogen_mass_fraction); + + /* Mean molecular mass for fully ionised gas */ + p->mu_ionised = 4. / (8. - 5. * (1. - p->hydrogen_mass_fraction)); + /* Read the artificial viscosity parameters from the file, if they exist */ p->viscosity.alpha = parser_get_opt_param_float( params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha); diff --git a/src/hydro_properties.h b/src/hydro_properties.h index b45b93192..5ee6a22d2 100644 --- a/src/hydro_properties.h +++ b/src/hydro_properties.h @@ -84,6 +84,12 @@ struct hydro_props { /*! Temperature of the neutral to ionized transition of Hydrogen */ float hydrogen_ionization_temperature; + /*! Mean molecular weight below hydrogen ionization temperature */ + float mu_neutral; + + /*! Mean molecular weight above hydrogen ionization temperature */ + float mu_ionised; + /*! Artificial viscosity parameters */ struct { /*! For the fixed, simple case. Also used to set the initial AV -- GitLab From 39d5daec43641842b691b5b89214b6bdc24d14eb Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Tue, 11 Dec 2018 10:10:02 +0000 Subject: [PATCH 17/22] Added calculation of temperature to the snapshots when using the constant Lambda cooling model. --- src/cooling/const_lambda/cooling.h | 43 +++++++++++++++++++++++++++ src/cooling/const_lambda/cooling_io.h | 18 +++++++++-- 2 files changed, 59 insertions(+), 2 deletions(-) diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h index 3c336060b..09c96413c 100644 --- a/src/cooling/const_lambda/cooling.h +++ b/src/cooling/const_lambda/cooling.h @@ -230,6 +230,49 @@ __attribute__((always_inline)) INLINE static void cooling_first_init_part( xp->cooling_data.radiated_energy = 0.f; } +/** + * @brief Compute the temperature of a #part based on the cooling function. + * + * @param phys_const #phys_const data structure. + * @param hydro_props The properties of the hydro scheme. + * @param us The internal system of units. + * @param cosmo #cosmology data structure. + * @param cooling #cooling_function_data struct. + * @param p #part data. + * @param xp Pointer to the #xpart data. + */ +INLINE static float cooling_get_temperature( + const struct phys_const* restrict phys_const, + const struct hydro_props* restrict hydro_props, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* restrict cooling, + const struct part* restrict p, const struct xpart* restrict xp) { + + /* Physical constants */ + const double m_H = phys_const->const_proton_mass; + const double k_B = phys_const->const_boltzmann_k; + + /* Gas properties */ + const double T_transition = hydro_props->hydrogen_ionization_temperature; + const double mu_neutral = hydro_props->mu_neutral; + const double mu_ionised = hydro_props->mu_ionised; + + /* Particle temperature */ + const double u = hydro_get_physical_internal_energy(p, xp, cosmo); + + /* Temperature over mean molecular weight */ + const double T_over_mu = hydro_gamma_minus_one * u * m_H / k_B; + + /* Are we above or below the HII -> HI transition? */ + if (T_over_mu > (T_transition + 1.) / mu_ionised) + return T_over_mu * mu_ionised; + else if (T_over_mu < (T_transition - 1.) / mu_neutral) + return T_over_mu * mu_neutral; + else + return T_transition; +} + /** * @brief Returns the total radiated energy by this particle. * diff --git a/src/cooling/const_lambda/cooling_io.h b/src/cooling/const_lambda/cooling_io.h index 0dca5011e..abd586b61 100644 --- a/src/cooling/const_lambda/cooling_io.h +++ b/src/cooling/const_lambda/cooling_io.h @@ -32,6 +32,7 @@ #include "../config.h" /* Local includes */ +#include "cooling.h" #include "io_properties.h" #ifdef HAVE_HDF5 @@ -49,11 +50,20 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( } #endif +INLINE static void convert_part_T(const struct engine* e, const struct part* p, + const struct xpart* xp, float* ret) { + + ret[0] = cooling_get_temperature(e->physical_constants, e->hydro_properties, + e->internal_units, e->cosmology, + e->cooling_func, p, xp); +} + /** * @brief Specifies which particle fields to write to a dataset * * Nothing to write for this scheme. * + * @param parts The particle array. * @param xparts The extended particle array. * @param list The list of i/o properties to write. * @param cooling The #cooling_function_data @@ -61,10 +71,14 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( * @return Returns the number of fields to write. */ __attribute__((always_inline)) INLINE static int cooling_write_particles( - const struct xpart* xparts, struct io_props* list, + const struct part* parts, const struct xpart* xparts, struct io_props* list, const struct cooling_function_data* cooling) { - return 0; + list[0] = io_make_output_field_convert_part("Temperature", FLOAT, 1, + UNIT_CONV_TEMPERATURE, parts, + xparts, convert_part_T); + + return 1; } #endif /* SWIFT_COOLING_CONST_LAMBDA_IO_H */ -- GitLab From 0cc541dc54bae69e5e9980708c0e6bebfd0dacc0 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Tue, 11 Dec 2018 10:18:56 +0000 Subject: [PATCH 18/22] Added calculation of temperature to the snapshots when using the constant du/dt cooling model. --- src/cooling/const_du/cooling.h | 43 +++++++++++++++++++++++++++++++ src/cooling/const_du/cooling_io.h | 19 ++++++++++++-- 2 files changed, 60 insertions(+), 2 deletions(-) diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h index dac92f098..8dc545f5b 100644 --- a/src/cooling/const_du/cooling.h +++ b/src/cooling/const_du/cooling.h @@ -164,6 +164,49 @@ __attribute__((always_inline)) INLINE static void cooling_first_init_part( xp->cooling_data.radiated_energy = 0.f; } +/** + * @brief Compute the temperature of a #part based on the cooling function. + * + * @param phys_const #phys_const data structure. + * @param hydro_props The properties of the hydro scheme. + * @param us The internal system of units. + * @param cosmo #cosmology data structure. + * @param cooling #cooling_function_data struct. + * @param p #part data. + * @param xp Pointer to the #xpart data. + */ +INLINE static float cooling_get_temperature( + const struct phys_const* restrict phys_const, + const struct hydro_props* restrict hydro_props, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* restrict cooling, + const struct part* restrict p, const struct xpart* restrict xp) { + + /* Physical constants */ + const double m_H = phys_const->const_proton_mass; + const double k_B = phys_const->const_boltzmann_k; + + /* Gas properties */ + const double T_transition = hydro_props->hydrogen_ionization_temperature; + const double mu_neutral = hydro_props->mu_neutral; + const double mu_ionised = hydro_props->mu_ionised; + + /* Particle temperature */ + const double u = hydro_get_physical_internal_energy(p, xp, cosmo); + + /* Temperature over mean molecular weight */ + const double T_over_mu = hydro_gamma_minus_one * u * m_H / k_B; + + /* Are we above or below the HII -> HI transition? */ + if (T_over_mu > (T_transition + 1.) / mu_ionised) + return T_over_mu * mu_ionised; + else if (T_over_mu < (T_transition - 1.) / mu_neutral) + return T_over_mu * mu_neutral; + else + return T_transition; +} + /** * @brief Returns the total radiated energy by this particle. * diff --git a/src/cooling/const_du/cooling_io.h b/src/cooling/const_du/cooling_io.h index f4a327f14..a60aa5d28 100644 --- a/src/cooling/const_du/cooling_io.h +++ b/src/cooling/const_du/cooling_io.h @@ -34,6 +34,7 @@ #include "../config.h" /* Local includes */ +#include "cooling.h" #include "io_properties.h" #ifdef HAVE_HDF5 @@ -50,9 +51,18 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( } #endif +INLINE static void convert_part_T(const struct engine* e, const struct part* p, + const struct xpart* xp, float* ret) { + + ret[0] = cooling_get_temperature(e->physical_constants, e->hydro_properties, + e->internal_units, e->cosmology, + e->cooling_func, p, xp); +} + /** * @brief Specifies which particle fields to write to a dataset * + * @param parts The particle array. * @param xparts The exended particle data array. * @param list The list of i/o properties to write. * @param cooling The #cooling_function_data @@ -60,9 +70,14 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( * @return Returns the number of fields to write. */ __attribute__((always_inline)) INLINE static int cooling_write_particles( - const struct xpart* xparts, struct io_props* list, + const struct part* parts, const struct xpart* xparts, struct io_props* list, const struct cooling_function_data* cooling) { - return 0; + + list[0] = io_make_output_field_convert_part("Temperature", FLOAT, 1, + UNIT_CONV_TEMPERATURE, parts, + xparts, convert_part_T); + + return 1; } #endif /* SWIFT_COOLING_CONST_DU_IO_H */ -- GitLab From 98629e45fb025c9ac65c6beb81b7102842fc7088 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Tue, 11 Dec 2018 10:28:13 +0000 Subject: [PATCH 19/22] Added calculation of temperature to the snapshots when using the no cooling model. --- src/cooling/none/cooling.h | 43 +++++++++++++++++++++++++++++++++++ src/cooling/none/cooling_io.h | 18 +++++++++++++-- 2 files changed, 59 insertions(+), 2 deletions(-) diff --git a/src/cooling/none/cooling.h b/src/cooling/none/cooling.h index 868bfad7f..579cf2ae9 100644 --- a/src/cooling/none/cooling.h +++ b/src/cooling/none/cooling.h @@ -119,6 +119,49 @@ __attribute__((always_inline)) INLINE static void cooling_first_init_part( const struct cooling_function_data* data, const struct part* restrict p, struct xpart* restrict xp) {} +/** + * @brief Compute the temperature of a #part based on the cooling function. + * + * @param phys_const #phys_const data structure. + * @param hydro_props The properties of the hydro scheme. + * @param us The internal system of units. + * @param cosmo #cosmology data structure. + * @param cooling #cooling_function_data struct. + * @param p #part data. + * @param xp Pointer to the #xpart data. + */ +INLINE static float cooling_get_temperature( + const struct phys_const* restrict phys_const, + const struct hydro_props* restrict hydro_props, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* restrict cooling, + const struct part* restrict p, const struct xpart* restrict xp) { + + /* Physical constants */ + const double m_H = phys_const->const_proton_mass; + const double k_B = phys_const->const_boltzmann_k; + + /* Gas properties */ + const double T_transition = hydro_props->hydrogen_ionization_temperature; + const double mu_neutral = hydro_props->mu_neutral; + const double mu_ionised = hydro_props->mu_ionised; + + /* Particle temperature */ + const double u = hydro_get_physical_internal_energy(p, xp, cosmo); + + /* Temperature over mean molecular weight */ + const double T_over_mu = hydro_gamma_minus_one * u * m_H / k_B; + + /* Are we above or below the HII -> HI transition? */ + if (T_over_mu > (T_transition + 1.) / mu_ionised) + return T_over_mu * mu_ionised; + else if (T_over_mu < (T_transition - 1.) / mu_neutral) + return T_over_mu * mu_neutral; + else + return T_transition; +} + /** * @brief Returns the total radiated energy by this particle. * diff --git a/src/cooling/none/cooling_io.h b/src/cooling/none/cooling_io.h index 518c16648..16b4b4ca2 100644 --- a/src/cooling/none/cooling_io.h +++ b/src/cooling/none/cooling_io.h @@ -23,6 +23,7 @@ #include "../config.h" /* Local includes */ +#include "cooling.h" #include "io_properties.h" #ifdef HAVE_HDF5 @@ -39,9 +40,18 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( } #endif +INLINE static void convert_part_T(const struct engine* e, const struct part* p, + const struct xpart* xp, float* ret) { + + ret[0] = cooling_get_temperature(e->physical_constants, e->hydro_properties, + e->internal_units, e->cosmology, + e->cooling_func, p, xp); +} + /** * @brief Specifies which particle fields to write to a dataset * + * @param parts The particle array. * @param xparts The extended particle array. * @param list The list of i/o properties to write. * @param cooling The #cooling_function_data @@ -49,9 +59,13 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( * @return Returns the number of fields to write. */ __attribute__((always_inline)) INLINE static int cooling_write_particles( - const struct xpart* xparts, struct io_props* list, + const struct part* parts, const struct xpart* xparts, struct io_props* list, const struct cooling_function_data* cooling) { - return 0; + + list[0] = io_make_output_field_convert_part("Temperature", FLOAT, 1, + UNIT_CONV_TEMPERATURE, parts, + xparts, convert_part_T); + return 1; } #endif /* SWIFT_COOLING_NONE_IO_H */ -- GitLab From a0a1eb90594fa8083ce9d64fbe398ed04dfd6252 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Tue, 11 Dec 2018 10:32:15 +0000 Subject: [PATCH 20/22] Also use the same function signature for the EAGLE cooling model when computing temperatures. --- src/cooling/EAGLE/cooling.c | 2 ++ src/cooling/EAGLE/cooling.h | 1 + src/cooling/EAGLE/cooling_io.h | 5 +++-- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c index cd3d80efb..d4755d055 100644 --- a/src/cooling/EAGLE/cooling.c +++ b/src/cooling/EAGLE/cooling.c @@ -644,6 +644,7 @@ __attribute__((always_inline)) INLINE void cooling_first_init_part( * The temperature returned is consistent with the cooling rates. * * @param phys_const #phys_const data structure. + * @param hydro_props The properties of the hydro scheme. * @param us The internal system of units. * @param cosmo #cosmology data structure. * @param cooling #cooling_function_data struct. @@ -652,6 +653,7 @@ __attribute__((always_inline)) INLINE void cooling_first_init_part( */ float cooling_get_temperature( const struct phys_const *restrict phys_const, + const struct hydro_props *restrict hydro_props, const struct unit_system *restrict us, const struct cosmology *restrict cosmo, const struct cooling_function_data *restrict cooling, diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h index 733749548..02ee37a41 100644 --- a/src/cooling/EAGLE/cooling.h +++ b/src/cooling/EAGLE/cooling.h @@ -61,6 +61,7 @@ void cooling_first_init_part( float cooling_get_temperature( const struct phys_const *restrict phys_const, + const struct hydro_props *restrict hydro_props, const struct unit_system *restrict us, const struct cosmology *restrict cosmo, const struct cooling_function_data *restrict cooling, diff --git a/src/cooling/EAGLE/cooling_io.h b/src/cooling/EAGLE/cooling_io.h index 86b7067bd..5508153af 100644 --- a/src/cooling/EAGLE/cooling_io.h +++ b/src/cooling/EAGLE/cooling_io.h @@ -44,8 +44,9 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( INLINE static void convert_part_T(const struct engine* e, const struct part* p, const struct xpart* xp, float* ret) { - ret[0] = cooling_get_temperature(e->physical_constants, e->internal_units, - e->cosmology, e->cooling_func, p, xp); + ret[0] = cooling_get_temperature(e->physical_constants, e->hydro_properties, + e->internal_units, e->cosmology, + e->cooling_func, p, xp); } /** -- GitLab From c3821386b7d70c0af470f346eb3799c6caa3da4c Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 12 Dec 2018 08:52:27 +0000 Subject: [PATCH 21/22] Added temperature snapshot output to the Compton cooling model. --- src/cooling/Compton/cooling.h | 43 +++++++++++++++++++++++++++ src/cooling/Compton/cooling_io.h | 20 ++++++++++--- src/cooling/const_lambda/cooling_io.h | 2 -- 3 files changed, 59 insertions(+), 6 deletions(-) diff --git a/src/cooling/Compton/cooling.h b/src/cooling/Compton/cooling.h index f440cd034..d6ec24f02 100644 --- a/src/cooling/Compton/cooling.h +++ b/src/cooling/Compton/cooling.h @@ -273,6 +273,49 @@ __attribute__((always_inline)) INLINE static void cooling_first_init_part( xp->cooling_data.radiated_energy = 0.f; } +/** + * @brief Compute the temperature of a #part based on the cooling function. + * + * @param phys_const #phys_const data structure. + * @param hydro_props The properties of the hydro scheme. + * @param us The internal system of units. + * @param cosmo #cosmology data structure. + * @param cooling #cooling_function_data struct. + * @param p #part data. + * @param xp Pointer to the #xpart data. + */ +INLINE static float cooling_get_temperature( + const struct phys_const* restrict phys_const, + const struct hydro_props* restrict hydro_props, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* restrict cooling, + const struct part* restrict p, const struct xpart* restrict xp) { + + /* Physical constants */ + const double m_H = phys_const->const_proton_mass; + const double k_B = phys_const->const_boltzmann_k; + + /* Gas properties */ + const double T_transition = hydro_props->hydrogen_ionization_temperature; + const double mu_neutral = hydro_props->mu_neutral; + const double mu_ionised = hydro_props->mu_ionised; + + /* Particle temperature */ + const double u = hydro_get_physical_internal_energy(p, xp, cosmo); + + /* Temperature over mean molecular weight */ + const double T_over_mu = hydro_gamma_minus_one * u * m_H / k_B; + + /* Are we above or below the HII -> HI transition? */ + if (T_over_mu > (T_transition + 1.) / mu_ionised) + return T_over_mu * mu_ionised; + else if (T_over_mu < (T_transition - 1.) / mu_neutral) + return T_over_mu * mu_neutral; + else + return T_transition; +} + /** * @brief Returns the total radiated energy by this particle. * diff --git a/src/cooling/Compton/cooling_io.h b/src/cooling/Compton/cooling_io.h index d020587c9..8fa3944ff 100644 --- a/src/cooling/Compton/cooling_io.h +++ b/src/cooling/Compton/cooling_io.h @@ -23,6 +23,7 @@ #include "../config.h" /* Local includes */ +#include "cooling.h" #include "io_properties.h" #ifdef HAVE_HDF5 @@ -41,11 +42,18 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( } #endif +INLINE static void convert_part_T(const struct engine* e, const struct part* p, + const struct xpart* xp, float* ret) { + + ret[0] = cooling_get_temperature(e->physical_constants, e->hydro_properties, + e->internal_units, e->cosmology, + e->cooling_func, p, xp); +} + /** * @brief Specifies which particle fields to write to a dataset * - * Nothing to write for this scheme. - * + * @param parts The particle array. * @param xparts The extended particle array. * @param list The list of i/o properties to write. * @param cooling The #cooling_function_data @@ -53,10 +61,14 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( * @return Returns the number of fields to write. */ __attribute__((always_inline)) INLINE static int cooling_write_particles( - const struct xpart* xparts, struct io_props* list, + const struct part* parts, const struct xpart* xparts, struct io_props* list, const struct cooling_function_data* cooling) { - return 0; + list[0] = io_make_output_field_convert_part("Temperature", FLOAT, 1, + UNIT_CONV_TEMPERATURE, parts, + xparts, convert_part_T); + + return 1; } #endif /* SWIFT_COOLING_IO_COMPTON_H */ diff --git a/src/cooling/const_lambda/cooling_io.h b/src/cooling/const_lambda/cooling_io.h index abd586b61..9437f0f94 100644 --- a/src/cooling/const_lambda/cooling_io.h +++ b/src/cooling/const_lambda/cooling_io.h @@ -61,8 +61,6 @@ INLINE static void convert_part_T(const struct engine* e, const struct part* p, /** * @brief Specifies which particle fields to write to a dataset * - * Nothing to write for this scheme. - * * @param parts The particle array. * @param xparts The extended particle array. * @param list The list of i/o properties to write. -- GitLab From e2f4cac4dfaed51ad46cab3a302d7510452c3cb6 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 12 Dec 2018 09:00:12 +0000 Subject: [PATCH 22/22] Updated the Compton cooling model to use the newly computed global quantities. --- src/cooling/Compton/cooling.h | 83 ++++++++++------------------------- 1 file changed, 23 insertions(+), 60 deletions(-) diff --git a/src/cooling/Compton/cooling.h b/src/cooling/Compton/cooling.h index d6ec24f02..77252140d 100644 --- a/src/cooling/Compton/cooling.h +++ b/src/cooling/Compton/cooling.h @@ -49,61 +49,11 @@ INLINE static void cooling_update(const struct cosmology* cosmo, // Add content if required. } -/** - * @brief Compute the mean molecular weight as a function of temperature for - * primordial gas. - * - * @param T The temperature of the gas [K]. - * @param H_mass_fraction The hydrogen mass fraction of the gas. - * @param T_transition The temperature of the transition from HII to HI [K]. - */ -__attribute__((always_inline, const)) INLINE static double -mean_molecular_weight(const double T, const double H_mass_fraction, - const double T_transition) { - - if (T > T_transition) - return 4. / (8. - 5. * (1. - H_mass_fraction)); - else - return 4. / (1. + 3. * H_mass_fraction); -} - -/** - * @brief Compute the temperature for a given internal energy per unit mass - * assuming primordial gas. - * - * @param u_cgs The internal energy per unit mass of the gas [erg * g^-1]. - * @param H_mass_fraction The hydrogen mass fraction of the gas. - * @param T_transition The temperature of the transition from HII to HI [K]. - * @param m_H_cgs The mass of the Hydorgen atom [g]. - * @param k_B_cgs The Boltzmann constant in cgs units [erg * K^-1] - * @return The temperature of the gas [K] - */ -__attribute__((always_inline, const)) INLINE static double -temperature_from_internal_energy(const double u_cgs, - const double H_mass_fraction, - const double T_transition, - const double m_H_cgs, const double k_B_cgs) { - - const double T_over_mu = hydro_gamma_minus_one * u_cgs * m_H_cgs / k_B_cgs; - - const double mu_high = - mean_molecular_weight(T_transition + 1., H_mass_fraction, T_transition); - const double mu_low = - mean_molecular_weight(T_transition - 1., H_mass_fraction, T_transition); - - if (T_over_mu > (T_transition + 1.) / mu_high) - return T_over_mu * mu_high; - else if (T_over_mu < (T_transition - 1.) / mu_low) - return T_over_mu * mu_low; - else - return T_transition; -} - /** * @brief Calculates du/dt in CGS units for a particle. * - * * @param cosmo The current cosmological model. + * @param phys_const The physical constants in internal units. * @param hydro_props The properties of the hydro scheme. * @param cooling The #cooling_function_data used in the run. * @param z The current redshift. @@ -113,7 +63,8 @@ temperature_from_internal_energy(const double u_cgs, * in cgs units [erg * g^-1 * s^-1]. */ __attribute__((always_inline)) INLINE static double Compton_cooling_rate_cgs( - const struct cosmology* cosmo, const struct hydro_props* hydro_props, + const struct cosmology* cosmo, const struct phys_const* restrict phys_const, + const struct hydro_props* hydro_props, const struct cooling_function_data* cooling, const double z, const double u, const struct part* p) { @@ -129,15 +80,27 @@ __attribute__((always_inline)) INLINE static double Compton_cooling_rate_cgs( /* CMB temperature at this redshift */ const double T_CMB = cooling->const_T_CMB_0 * zp1; + /* Physical constants */ + const double m_H = phys_const->const_proton_mass; + const double k_B = phys_const->const_boltzmann_k; + /* Gas properties */ - const double H_mass_fraction = hydro_props->hydrogen_mass_fraction; const double T_transition = hydro_props->hydrogen_ionization_temperature; + const double mu_neutral = hydro_props->mu_neutral; + const double mu_ionised = hydro_props->mu_ionised; - /* Particle temperature */ - const double u_cgs = u * cooling->conv_factor_energy_to_cgs; - const double T = temperature_from_internal_energy(u_cgs, H_mass_fraction, - T_transition, 1., 1.); - // MATTHIEU: to do: get H mass in cgs and k_B in cgs. + /* Temperature over mean molecular weight */ + const double T_over_mu = hydro_gamma_minus_one * u * m_H / k_B; + + double T; + + /* Are we above or below the HII -> HI transition? */ + if (T_over_mu > (T_transition + 1.) / mu_ionised) + T = T_over_mu * mu_ionised; + else if (T_over_mu < (T_transition - 1.) / mu_neutral) + T = T_over_mu * mu_neutral; + else + T = T_transition; /* Electron abundance */ double electron_abundance = 0.; // MATTHIEU: To do: compute X_e @@ -189,8 +152,8 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( const float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo); /* Calculate cooling du_dt (in cgs units) */ - const double cooling_du_dt_cgs = - Compton_cooling_rate_cgs(cosmo, hydro_props, cooling, cosmo->z, u_old, p); + const double cooling_du_dt_cgs = Compton_cooling_rate_cgs( + cosmo, phys_const, hydro_props, cooling, cosmo->z, u_old, p); /* Convert to internal units */ float cooling_du_dt = -- GitLab