From bc6998e21bc7dc10eece24b6df5db0025d05906d Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Mon, 8 Oct 2018 22:11:43 +0200
Subject: [PATCH] Added function to compute T from u in the Compton cooling
 example.

---
 .../plotTempEvolution.py                      |  2 +-
 src/cooling/Compton/cooling.h                 | 67 +++++++++++++++++--
 src/cooling/Compton/cooling_struct.h          |  3 +
 3 files changed, 67 insertions(+), 5 deletions(-)

diff --git a/examples/SmallCosmoVolume_cooling/plotTempEvolution.py b/examples/SmallCosmoVolume_cooling/plotTempEvolution.py
index 790e80bf9b..d3e4e9047e 100644
--- a/examples/SmallCosmoVolume_cooling/plotTempEvolution.py
+++ b/examples/SmallCosmoVolume_cooling/plotTempEvolution.py
@@ -85,7 +85,7 @@ unit_length_in_si = 0.01 * unit_length_in_cgs
 unit_mass_in_si = 0.001 * unit_mass_in_cgs
 unit_time_in_si = unit_time_in_cgs
 
-# Primoridal ean molecular weight as a function of temperature
+# Primoridal mean molecular weight as a function of temperature
 def mu(T, H_frac=H_mass_fraction, T_trans=H_transition_temp):
     if T > T_trans:
         return 4. / (8. - 5. * (1. - H_frac))
diff --git a/src/cooling/Compton/cooling.h b/src/cooling/Compton/cooling.h
index 03cd21a4f3..c8da002655 100644
--- a/src/cooling/Compton/cooling.h
+++ b/src/cooling/Compton/cooling.h
@@ -35,6 +35,56 @@
 #include "physical_constants.h"
 #include "units.h"
 
+/**
+ * @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_trans 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_trans 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.
  *
@@ -42,14 +92,15 @@
  * @param cosmo The current cosmological model.
  * @param hydro_props The properties of the hydro scheme.
  * @param cooling The #cooling_function_data used in the run.
- * @param z The current redshift
+ * @param z The current redshift.
+ * @param u The current internal energy in internal units.
  * @param p Pointer to the particle data.
  * @return The change in energy per unit mass due to cooling for this particle
  * 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 cooling_function_data* cooling, const double z,
+    const struct cooling_function_data* cooling, const double z, const double u,
     const struct part* p) {
 
   /* Get particle density */
@@ -64,8 +115,14 @@ __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;
 
+  /* Gas properties */
+  const double H_mass_fraction = hydro_props->hydrogen_mass_fraction;
+  const double T_transition = hydro_props->hydrogen_ionization_temperature;
+
   /* Particle temperature */
-  const double T = 1;
+  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
 
   /* Temperature difference with the CMB */
   const double delta_T = T - T_CMB;
@@ -115,7 +172,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
 
   /* Calculate cooling du_dt (in cgs units) */
   const double cooling_du_dt_cgs =
-      Compton_cooling_rate_cgs(cosmo, hydro_props, cooling, cosmo->z, p);
+      Compton_cooling_rate_cgs(cosmo, hydro_props, cooling, cosmo->z, u_old, p);
 
   /* Convert to internal units */
   float cooling_du_dt =
@@ -225,6 +282,8 @@ static INLINE void cooling_init_backend(struct swift_params* parameter_file,
   /* Some useful conversion values */
   cooling->conv_factor_density_to_cgs =
       units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+  cooling->conv_factor_energy_to_cgs =
+      units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
   cooling->conv_factor_energy_rate_from_cgs =
       units_cgs_conversion_factor(us, UNIT_CONV_TIME) /
       units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
diff --git a/src/cooling/Compton/cooling_struct.h b/src/cooling/Compton/cooling_struct.h
index 16c7815818..7f1cb4c44b 100644
--- a/src/cooling/Compton/cooling_struct.h
+++ b/src/cooling/Compton/cooling_struct.h
@@ -36,6 +36,9 @@ struct cooling_function_data {
   /*! Conversion factor from internal units to cgs for density */
   double conv_factor_density_to_cgs;
 
+  /*! Conversion factor from internal units to cgs for internal energy */
+  double conv_factor_energy_to_cgs;
+
   /*! Conversion factor from internal units from cgs for internal energy
    * derivative */
   double conv_factor_energy_rate_from_cgs;
-- 
GitLab