Include chemistry in sph

src/cooling/grackle/cooling_struct.h

@@ -52,7 +52,9 @@ struct cooling_function_data {
 };
 
 /**
- * @brief Properties of the cooling stored in the particle data
+ * @brief Properties of the cooling stored in the extra particle data
+ *
+ * theses data are not processed during the SPH density/force.
  */
 struct cooling_xpart_data {
 
@@ -60,4 +62,13 @@ struct cooling_xpart_data {
   float radiated_energy;
 };
 
+/**
+ * @brief Properties of the cooling stored in the particle data
+ */
+struct cooling_part_data {
+
+  /*! Quick example */
+  float He_density;
+};
+
 #endif /* SWIFT_COOLING_STRUCT_NONE_H */

src/hydro/Gadget2/hydro_iact.h

@@ -32,6 +32,7 @@
  * Gadget-2 tree-code neighbours search.
  */
 
+#include "cooling.h"
 #include "cache.h"
 #include "minmax.h"
 
@@ -104,6 +105,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   pj->density.rot_v[1] += facj * curlvr[1];
   pj->density.rot_v[2] += facj * curlvr[2];
 
+  cooling_density_iact(wi, wj, pi, pj);
+
 #ifdef DEBUG_INTERACTIONS_SPH
   /* Update ngb counters */
   if (pi->num_ngb_density < MAX_NUM_OF_NEIGHBOURS)

src/hydro/Gadget2/hydro_io.h

@@ -23,6 +23,7 @@
 #include "hydro.h"
 #include "io_properties.h"
 #include "kernel_hydro.h"
+#include "cooling_io.h"
 
 /**
  * @brief Specifies which particle fields to read from a dataset
@@ -53,6 +54,8 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
                                 UNIT_CONV_ACCELERATION, parts, a_hydro);
   list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
                                 UNIT_CONV_DENSITY, parts, rho);
+
+  cooling_read_particles(parts, list, num_fields);
 }
 
 void convert_u(const struct engine* e, const struct part* p, float* ret) {
@@ -117,7 +120,9 @@ void hydro_write_particles(const struct part* parts, struct io_props* list,
                                               parts, convert_u);
   list[9] = io_make_output_field_convert_part(
       "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
+
 #ifdef DEBUG_INTERACTIONS_SPH
+  
   list[10] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS,
                                   parts, num_ngb_density);
   list[11] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS,
@@ -129,6 +134,8 @@ void hydro_write_particles(const struct part* parts, struct io_props* list,
       io_make_output_field("Ids_ngb_force", LONGLONG, MAX_NUM_OF_NEIGHBOURS,
                            UNIT_CONV_NO_UNITS, parts, ids_ngbs_force);
 #endif
+
+  cooling_write_particles(parts, list, num_fields);
 }
 
 /**
@@ -145,6 +152,8 @@ void writeSPHflavour(hid_t h_grpsph) {
       "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
   io_write_attribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
   io_write_attribute_f(h_grpsph, "Viscosity beta", 3.f);
+
+  writeCoolingFlavor(h_grpsph);
 }
 
 /**

src/hydro/Gadget2/hydro_part.h

@@ -133,6 +133,8 @@ struct part {
   /* Time-step length */
   timebin_t time_bin;
 
+  struct cooling_part_data cooling_data;
+  
 #ifdef SWIFT_DEBUG_CHECKS
 
   /* Time of the last drift */