From 74e7fc6d8e17e71967c3ed6b74d72d0ee4d17c10 Mon Sep 17 00:00:00 2001
From: lhausamm <loic_hausammann@hotmail.com>
Date: Mon, 8 Jan 2018 15:55:39 +0100
Subject: [PATCH] Include chemistry in sph

---
 src/cooling/grackle/cooling_struct.h | 13 ++++++++++++-
 src/hydro/Gadget2/hydro_iact.h       |  3 +++
 src/hydro/Gadget2/hydro_io.h         |  9 +++++++++
 src/hydro/Gadget2/hydro_part.h       |  2 ++
 4 files changed, 26 insertions(+), 1 deletion(-)

diff --git a/src/cooling/grackle/cooling_struct.h b/src/cooling/grackle/cooling_struct.h
index ddc07ed70e..ec1bd6cb8f 100644
--- a/src/cooling/grackle/cooling_struct.h
+++ b/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 */
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index d30fd7b803..29a86da3f9 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/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)
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 34a9f33daf..7be629f488 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/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);
 }
 
 /**
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index c7deaf7360..4bfdbcf7ff 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/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 */
-- 
GitLab