Implement a way to set an initial temperature to the gas when specified in ICs.

examples/main.c

@@ -609,7 +609,8 @@ int main(int argc, char *argv[]) {
     if (myrank == 0 && with_cosmology) cosmology_print(&cosmo);
 
     /* Initialise the hydro properties */
-    if (with_hydro) hydro_props_init(&hydro_properties, params);
+    if (with_hydro)
+      hydro_props_init(&hydro_properties, &prog_const, &us, params);
     if (with_hydro) eos_init(&eos, params);
 
     /* Initialise the gravity properties */

examples/parameter_example.yml

@@ -30,7 +30,9 @@ SPH:
   h_max:                 10.      # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
   max_volume_change:     1.4      # (Optional) Maximal allowed change of kernel volume over one time-step.
   max_ghost_iterations:  30       # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
-
+  initial_temperature:   100      # (Optional) Initial temperature (in internal units) to set the gass particles to. Value is ignored if set to 0.
+  H_mass_fraction:       0.76     # (Optional) Hydrogen mass fraction used for initial conversion from temp to internal energy.
+  
 # Parameters for the self-gravity scheme
 Gravity:
   eta:          0.025    # Constant dimensionless multiplier for time integration.

src/hydro/Gadget2/hydro.h

@@ -598,4 +598,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
   hydro_init_part(p, NULL);
 }
 
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+  p->entropy = u_init;
+}
+
 #endif /* SWIFT_GADGET2_HYDRO_H */

src/hydro_properties.c

@@ -36,8 +36,12 @@
 #define hydro_props_default_volume_change 1.4f
 #define hydro_props_default_h_max FLT_MAX
 #define hydro_props_default_h_tolerance 1e-4
+#define hydro_props_default_init_temp 0.f
+#define hydro_props_default_H_fraction 0.76
 
 void hydro_props_init(struct hydro_props *p,
+                      const struct phys_const *phys_const,
+                      const struct unit_system *us,
                       const struct swift_params *params) {
 
   /* Kernel properties */
@@ -73,6 +77,30 @@ void hydro_props_init(struct hydro_props *p,
   const float max_volume_change = parser_get_opt_param_float(
       params, "SPH:max_volume_change", hydro_props_default_volume_change);
   p->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
+
+  /* Initial temperature */
+  p->initial_temperature = parser_get_opt_param_float(
+      params, "SPH:initial_temperature", hydro_props_default_init_temp);
+
+  /* Hydrogen mass fraction */
+  p->hydrogen_mass_fraction = parser_get_opt_param_double(
+      params, "SPH:H_mass_fraction", hydro_props_default_H_fraction);
+
+  /* Compute the initial energy (Note the temp. read is in internal units) */
+  double u_init = phys_const->const_boltzmann_k / phys_const->const_proton_mass;
+  u_init *= p->initial_temperature;
+  u_init *= hydro_one_over_gamma_minus_one;
+
+  /* Correct for hydrogen mass fraction */
+  double mean_molecular_weight;
+  if (p->initial_temperature *
+          units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE) >
+      1e4)
+    mean_molecular_weight = 4. / (8. - 5. * (1. - p->hydrogen_mass_fraction));
+  else
+    mean_molecular_weight = 4. / (1. + 3. * p->hydrogen_mass_fraction);
+
+  p->initial_internal_energy = u_init / mean_molecular_weight;
 }
 
 void hydro_props_print(const struct hydro_props *p) {
@@ -103,6 +131,9 @@ void hydro_props_print(const struct hydro_props *p) {
   if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
     message("Maximal iterations in ghost task set to %d (default is %d)",
             p->max_smoothing_iterations, hydro_props_default_max_iterations);
+
+  if (p->initial_temperature != hydro_props_default_init_temp)
+    message("Initial gas temperature set to %f", p->initial_temperature);
 }
 
 #if defined(HAVE_HDF5)
@@ -125,6 +156,11 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
                        pow_dimension(expf(p->log_max_h_change)));
   io_write_attribute_i(h_grpsph, "Max ghost iterations",
                        p->max_smoothing_iterations);
+  io_write_attribute_f(h_grpsph, "Initial temperature", p->initial_temperature);
+  io_write_attribute_f(h_grpsph, "Initial energy per unit mass",
+                       p->initial_internal_energy);
+  io_write_attribute_f(h_grpsph, "Hydrogen mass fraction",
+                       p->hydrogen_mass_fraction);
 }
 #endif
 

src/hydro_properties.h

@@ -33,7 +33,9 @@
 
 /* Local includes. */
 #include "parser.h"
+#include "physical_constants.h"
 #include "restart.h"
+#include "units.h"
 
 /**
  * @brief Contains all the constants and parameters of the hydro scheme
@@ -63,10 +65,22 @@ struct hydro_props {
 
   /*! Maximal change of h over one time-step */
   float log_max_h_change;
+
+  /*! Initial temperature */
+  float initial_temperature;
+
+  /*! Initial internal energy per unit mass */
+  float initial_internal_energy;
+
+  /*! Primoridal hydrogen mass fraction for initial energy conversion */
+  float hydrogen_mass_fraction;
 };
 
 void hydro_props_print(const struct hydro_props *p);
-void hydro_props_init(struct hydro_props *p, const struct swift_params *params);
+void hydro_props_init(struct hydro_props *p,
+                      const struct phys_const *phys_const,
+                      const struct unit_system *us,
+                      const struct swift_params *params);
 
 #if defined(HAVE_HDF5)
 void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p);

src/space.c

@@ -2646,9 +2646,13 @@ void space_first_init_parts(struct space *s,
   const size_t nr_parts = s->nr_parts;
   struct part *restrict p = s->parts;
   struct xpart *restrict xp = s->xparts;
+
   const struct cosmology *cosmo = s->e->cosmology;
   const float a_factor_vel = cosmo->a * cosmo->a;
 
+  const struct hydro_props *hydro_props = s->e->hydro_properties;
+  const float u_init = hydro_props->initial_internal_energy;
+
   for (size_t i = 0; i < nr_parts; ++i) {
 
     /* Convert velocities to internal units */
@@ -2668,6 +2672,9 @@ void space_first_init_parts(struct space *s,
 
     hydro_first_init_part(&p[i], &xp[i]);
 
+    /* Overwrite the internal energy? */
+    if (u_init > 0.f) hydro_set_init_internal_energy(&p[i], u_init);
+
     /* Also initialise the chemistry */
     chemistry_first_init_part(&p[i], &xp[i], chemistry);
 
@@ -3254,6 +3261,11 @@ void space_generate_gas(struct space *s, const struct cosmology *cosmo,
     /* Set the smoothing length to the mean inter-particle separation */
     p->h = d;
   }
+
+  /* Replace the content of the space */
+  free(s->gparts);
+  s->parts = parts;
+  s->gparts = gparts;
 }
 
 /**