Grackle chemistry --> GEAR chemistry. Moved the output of subgrid model to a...

Grackle chemistry --> GEAR chemistry. Moved the output of subgrid model to a separate part of the HDF5 file.

configure.ac

@@ -935,7 +935,7 @@ esac
 #  chemistry function
 AC_ARG_WITH([chemistry],
    [AS_HELP_STRING([--with-chemistry=<function>],
-      [chemistry function @<:@none, grackle default: none@:>@]
+      [chemistry function @<:@none, gear default: none@:>@]
    )],
    [with_chemistry="$withval"],
    [with_chemistry="none"]
@@ -944,8 +944,8 @@ case "$with_chemistry" in
    none)
       AC_DEFINE([CHEMISTRY_NONE], [1], [No chemistry function])
    ;; 
-   grackle)
-      AC_DEFINE([CHEMISTRY_GRACKLE], [1], [Chemistry compatible with the grackle cooling])
+   gear)
+      AC_DEFINE([CHEMISTRY_GEAR], [1], [Chemistry taken from the GEAR model])
    ;; 
    *)
       AC_MSG_ERROR([Unknown chemistry function: $with_chemistry])

doc/Doxyfile.in

@@ -765,7 +765,8 @@ INPUT		       += @top_srcdir@/src/gravity/Default
 INPUT		       += @top_srcdir@/src/stars/Default
 INPUT		       += @top_srcdir@/src/riemann
 INPUT		       += @top_srcdir@/src/potential/point_mass
-INPUT		       += @top_srcdir@/src/cooling/const_du
+INPUT		       += @top_srcdir@/src/cooling/none
+INPUT		       += @top_srcdir@/src/chemistry/none
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses

examples/main.c

@@ -612,6 +612,9 @@ int main(int argc, char *argv[]) {
   if (with_cooling) cooling_init(params, &us, &prog_const, &cooling_func);
   if (with_cooling && myrank == 0) cooling_print(&cooling_func);
 
+  /* Initialise the chemistry */
+  chemistry_print();
+
   /* Initialise the feedback properties */
   struct sourceterms sourceterms;
   if (with_sourceterms) sourceterms_init(params, &us, &sourceterms);

src/Makefile.am

@@ -123,9 +123,9 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  chemistry/none/chemistry.h \
 		 chemistry/none/chemistry_io.h \
 		 chemistry/none/chemistry_struct.h \
-                 chemistry/grackle/chemistry.h \
-		 chemistry/grackle/chemistry_io.h \
-		 chemistry/grackle/chemistry_struct.h
+                 chemistry/gear/chemistry.h \
+		 chemistry/gear/chemistry_io.h \
+		 chemistry/gear/chemistry_struct.h
 
 
 # Sources and flags for regular library

src/chemistry.c

@@ -31,24 +31,17 @@
  * @param parameter_file The parsed parameter file.
  * @param us The current internal system of units.
  * @param phys_const The physical constants in internal units.
- * @param chemistry The chemistry properties to initialize
  */
 void chemistry_init(const struct swift_params* parameter_file,
                     const struct unit_system* us,
-                    const struct phys_const* phys_const,
-                    struct chemistry_part_data* chemistry) {
+                    const struct phys_const* phys_const) {
 
-  chemistry_init_backend(parameter_file, us, phys_const, chemistry);
+  chemistry_init_backend(parameter_file, us, phys_const);
 }
 
 /**
  * @brief Prints the properties of the chemistry model to stdout.
  *
- * Calls chemistry_print_backend for the chosen chemistry function.
- *
- * @param chemistry The properties of the chemistry function.
+ * Calls chemistry_print_backend for the chosen chemistry model.
  */
-void chemistry_print(const struct chemistry_part_data* chemistry) {
-
-  chemistry_print_backend(chemistry);
-}
+void chemistry_print() { chemistry_print_backend(); }

src/chemistry.h

@@ -31,8 +31,8 @@
 /* Import the right chemistry definition */
 #if defined(CHEMISTRY_NONE)
 #include "./chemistry/none/chemistry.h"
-#elif defined(CHEMISTRY_GRACKLE)
-#include "./chemistry/grackle/chemistry.h"
+#elif defined(CHEMISTRY_GEAR)
+#include "./chemistry/gear/chemistry.h"
 #else
 #error "Invalid choice of chemistry function."
 #endif
@@ -40,9 +40,8 @@
 /* Common functions */
 void chemistry_init(const struct swift_params* parameter_file,
                     const struct unit_system* us,
-                    const struct phys_const* phys_const,
-                    struct chemistry_part_data* chem);
+                    const struct phys_const* phys_const);
 
-void chemistry_print(const struct chemistry_part_data* chem);
+void chemistry_print();
 
 #endif /* SWIFT_CHEMISTRY_H */

src/chemistry/grackle/chemistry.h → src/chemistry/gear/chemistry.h

@@ -16,11 +16,11 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef SWIFT_CHEMISTRY_NONE_H
-#define SWIFT_CHEMISTRY_NONE_H
+#ifndef SWIFT_CHEMISTRY_GEAR_H
+#define SWIFT_CHEMISTRY_GEAR_H
 
 /**
- * @file src/chemistry/none/chemistry.h
+ * @file src/chemistry/gear/chemistry.h
  * @brief Empty infrastructure for the cases without chemistry function
  */
 
@@ -57,22 +57,17 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part(
  * @param parameter_file The parsed parameter file.
  * @param us The current internal system of units.
  * @param phys_const The physical constants in internal units.
- * @param chemistry The chemistry properties to initialize
  */
 static INLINE void chemistry_init_backend(
     const struct swift_params* parameter_file, const struct unit_system* us,
-    const struct phys_const* phys_const,
-    struct chemistry_part_data* chemistry) {}
+    const struct phys_const* phys_const) {}
 
 /**
  * @brief Prints the properties of the chemistry model to stdout.
- *
- * @param chemistry The properties of the chemistry function.
  */
-static INLINE void chemistry_print_backend(
-    const struct chemistry_part_data* chemistry) {
+static INLINE void chemistry_print_backend() {
 
-  message("Chemistry function is 'grackle'.");
+  message("Chemistry function is 'gear'.");
 }
 
-#endif /* SWIFT_CHEMISTRY_NONE_H */
+#endif /* SWIFT_CHEMISTRY_GEAR_H */

src/chemistry/grackle/chemistry_io.h → src/chemistry/gear/chemistry_io.h

@@ -16,6 +16,8 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
+#ifndef SWIFT_CHEMISTRY_IO_GEAR_H
+#define SWIFT_CHEMISTRY_IO_GEAR_H
 
 #include "io_properties.h"
 
@@ -24,18 +26,17 @@
  *
  * @param parts The particle array.
  * @param list The list of i/o properties to read.
- * @param num_fields The number of i/o fields to read.
+ *
+ * @return Returns the number of fields to read.
  */
-void chemistry_read_particles(struct part* parts, struct io_props* list,
-                              int* num_fields) {
-
-  list += *num_fields;
-  *num_fields += 1;
+int chemistry_read_particles(struct part* parts, struct io_props* list) {
 
   /* List what we want to read */
   list[0] =
       io_make_input_field("HeDensity", FLOAT, 1, COMPULSORY, UNIT_CONV_DENSITY,
                           parts, chemistry_data.he_density);
+
+  return 1;
 }
 
 /**
@@ -43,16 +44,16 @@ void chemistry_read_particles(struct part* parts, struct io_props* list,
  *
  * @param parts The particle array.
  * @param list The list of i/o properties to write.
- * @param num_fields The number of i/o fields to write.
+ *
+ * @return Returns the number of fields to write.
  */
-void chemistry_write_particles(const struct part* parts, struct io_props* list,
-                               int* num_fields) {
-
-  list += *num_fields;
-  *num_fields += 1;
+int chemistry_write_particles(const struct part* parts, struct io_props* list) {
 
+  /* List what we want to write */
   list[0] = io_make_output_field("HeDensity", FLOAT, 1, UNIT_CONV_DENSITY,
                                  parts, chemistry_data.he_density);
+
+  return 1;
 }
 
 /**
@@ -61,5 +62,7 @@ void chemistry_write_particles(const struct part* parts, struct io_props* list,
  */
 void writeChemistryFlavor(hid_t h_grpsph) {
 
-  io_write_attribute_s(h_grpsph, "Chemistry Model", "Grackle");
+  io_write_attribute_s(h_grpsph, "Chemistry Model", "GEAR");
 }
+
+#endif /* SWIFT_CHEMISTRY_IO_GEAR_H */

src/chemistry/grackle/chemistry_struct.h → src/chemistry/gear/chemistry_struct.h

@@ -16,11 +16,11 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef SWIFT_CHEMISTRY_STRUCT_NONE_H
-#define SWIFT_CHEMISTRY_STRUCT_NONE_H
+#ifndef SWIFT_CHEMISTRY_STRUCT_GEAR_H
+#define SWIFT_CHEMISTRY_STRUCT_GEAR_H
 
 /**
- * @file src/chemistry/none/chemistry_struct.h
+ * @file src/chemistry/gear/chemistry_struct.h
  * @brief Empty infrastructure for the cases without chemistry function
  */
 
@@ -31,4 +31,4 @@ struct chemistry_part_data {
   float he_density;
 };
 
-#endif /* SWIFT_CHEMISTRY_STRUCT_NONE_H */
+#endif /* SWIFT_CHEMISTRY_STRUCT_GEAR_H */

src/chemistry/none/chemistry.h

@@ -57,20 +57,15 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part(
  * @param parameter_file The parsed parameter file.
  * @param us The current internal system of units.
  * @param phys_const The physical constants in internal units.
- * @param chemistry The chemistry properties to initialize
  */
 static INLINE void chemistry_init_backend(
     const struct swift_params* parameter_file, const struct unit_system* us,
-    const struct phys_const* phys_const,
-    struct chemistry_part_data* chemistry) {}
+    const struct phys_const* phys_const) {}
 
 /**
  * @brief Prints the properties of the chemistry model to stdout.
- *
- * @param chemistry The properties of the chemistry function.
  */
-static INLINE void chemistry_print_backend(
-    const struct chemistry_part_data* chemistry) {
+static INLINE void chemistry_print_backend() {
 
   message("Chemistry function is 'No chemistry'.");
 }

src/chemistry/none/chemistry_io.h

@@ -16,6 +16,8 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
+#ifndef SWIFT_CHEMISTRY_IO_NONE_H
+#define SWIFT_CHEMISTRY_IO_NONE_H
 
 #include "io_properties.h"
 
@@ -24,12 +26,15 @@
  *
  * @param parts The particle array.
  * @param list The list of i/o properties to read.
- * @param num_fields The number of i/o fields to read.
+ *
+ * @return Returns the number of fields to write.
  */
-void chemistry_read_particles(struct part* parts, struct io_props* list,
-                              int* num_fields) {
+int chemistry_read_particles(struct part* parts, struct io_props* list) {
+
+  /* update list according to hydro_io */
 
-  /* update num_fields and list according to hydro_io */
+  /* Return the number of fields to read */
+  return 0;
 }
 
 /**
@@ -37,12 +42,15 @@ void chemistry_read_particles(struct part* parts, struct io_props* list,
  *
  * @param parts The particle array.
  * @param list The list of i/o properties to write.
- * @param num_fields The number of i/o fields to write.
+ *
+ * @return Returns the number of fields to write.
  */
-void chemistry_write_particles(const struct part* parts, struct io_props* list,
-                               int* num_fields) {
+int chemistry_write_particles(const struct part* parts, struct io_props* list) {
 
-  /* update num_fields and list according to hydro_io */
+  /* update list according to hydro_io */
+
+  /* Return the number of fields to write */
+  return 0;
 }
 
 /**
@@ -53,3 +61,5 @@ void writeChemistryFlavor(hid_t h_grpsph) {
 
   io_write_attribute_s(h_grpsph, "Chemistry Model", "None");
 }
+
+#endif /* SWIFT_CHEMISTRY_IO_NONE_H */

src/chemistry_io.h

@@ -25,8 +25,8 @@
 /* Import the right functions */
 #if defined(CHEMISTRY_NONE)
 #include "./chemistry/none/chemistry_io.h"
-#elif defined(CHEMISTRY_GRACKLE)
-#include "./chemistry/grackle/chemistry_io.h"
+#elif defined(CHEMISTRY_GEAR)
+#include "./chemistry/gear/chemistry_io.h"
 #else
 #error "Invalid choice of chemistry function."
 #endif

src/chemistry_struct.h

@@ -30,8 +30,8 @@
 /* Import the right chemistry definition */
 #if defined(CHEMISTRY_NONE)
 #include "./chemistry/none/chemistry_struct.h"
-#elif defined(CHEMISTRY_GRACKLE)
-#include "./chemistry/grackle/chemistry_struct.h"
+#elif defined(CHEMISTRY_GEAR)
+#include "./chemistry/gear/chemistry_struct.h"
 #else
 #error "Invalid choice of chemistry function."
 #endif

src/cooling/const_du/cooling.h

@@ -51,7 +51,7 @@
 __attribute__((always_inline)) INLINE static void writeCoolingFlavor(
     hid_t h_grpsph) {
 
-  io_write_attribute_s(h_grpsph, "Cooling Model", "const_du");
+  io_write_attribute_s(h_grpsph, "Cooling Model", "Constant du/dt");
 }
 
 /**

src/cooling/const_lambda/cooling.h

@@ -44,7 +44,7 @@
 __attribute__((always_inline)) INLINE static void writeCoolingFlavor(
     hid_t h_grpsph) {
 
-  io_write_attribute_s(h_grpsph, "Cooling Model", "const_lambda");
+  io_write_attribute_s(h_grpsph, "Cooling Model", "Constant Lambda");
 }
 
 /**

src/cooling/grackle/cooling_struct.h

@@ -53,8 +53,6 @@ struct cooling_function_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 {
 

src/hydro/Gadget2/hydro_io.h

@@ -36,8 +36,6 @@
 void hydro_read_particles(struct part* parts, struct io_props* list,
                           int* num_fields) {
 
-  *num_fields = 8;
-
   /* List what we want to read */
   list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
                                 UNIT_CONV_LENGTH, parts, x);
@@ -56,7 +54,13 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
   list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
                                 UNIT_CONV_DENSITY, parts, rho);
 
-  chemistry_read_particles(parts, list, num_fields);
+  *num_fields = 8;
+  list += *num_fields;
+
+  /* Read in chemistry information */
+  const int num_chem_fields = chemistry_read_particles(parts, list);
+  *num_fields += num_chem_fields;
+  list += num_chem_fields;
 }
 
 void convert_u(const struct engine* e, const struct part* p, float* ret) {
@@ -93,12 +97,6 @@ void convert_part_pos(const struct engine* e, const struct part* p,
 void hydro_write_particles(const struct part* parts, struct io_props* list,
                            int* num_fields) {
 
-  *num_fields = 10;
-
-#ifdef DEBUG_INTERACTIONS_SPH
-  *num_fields += 4;
-#endif
-
   /* List what we want to write */
   list[0] = io_make_output_field_convert_part(
       "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
@@ -122,21 +120,30 @@ void hydro_write_particles(const struct part* parts, struct io_props* list,
   list[9] = io_make_output_field_convert_part(
       "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
 
+  *num_fields = 10;
+  list += *num_fields;
+
+  /* Write some chemistry information */
+  const int num_chem_fields = chemistry_write_particles(parts, list);
+  *num_fields += num_chem_fields;
+  list += num_chem_fields;
+
 #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,
-                                  parts, num_ngb_force);
-  list[12] =
+  list[0] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS,
+                                 parts, num_ngb_density);
+  list[1] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS,
+                                 parts, num_ngb_force);
+  list[2] =
       io_make_output_field("Ids_ngb_density", LONGLONG, MAX_NUM_OF_NEIGHBOURS,
                            UNIT_CONV_NO_UNITS, parts, ids_ngbs_density);
-  list[13] =
+  list[3] =
       io_make_output_field("Ids_ngb_force", LONGLONG, MAX_NUM_OF_NEIGHBOURS,
                            UNIT_CONV_NO_UNITS, parts, ids_ngbs_force);
-#endif
 
-  chemistry_write_particles(parts, list, num_fields);
+  *num_fields += 4;
+  list += 4;
+#endif
 }
 
 /**
@@ -153,9 +160,6 @@ 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);
-  writeChemistryFlavor(h_grpsph);
 }
 
 /**

src/hydro/Gadget2/hydro_part.h

@@ -131,11 +131,12 @@ struct part {
     } force;
   };
 
+  /* Chemistry information */
+  struct chemistry_part_data chemistry_data;
+
   /* Time-step length */
   timebin_t time_bin;
 
-  struct chemistry_part_data chemistry_data;
-
 #ifdef SWIFT_DEBUG_CHECKS
 
   /* Time of the last drift */

src/parallel_io.c

@@ -961,6 +961,14 @@ void write_output_parallel(struct engine* e, const char* baseName,
     H5Gclose(h_grp);
   }
 
+  /* Print the subgrid parameters */
+  h_grp = H5Gcreate(h_file, "/SubgridScheme", H5P_DEFAULT, H5P_DEFAULT,
+                    H5P_DEFAULT);
+  if (h_grp < 0) error("Error while creating subgrid group");
+  writeCoolingFlavor(h_grp);
+  writeChemistryFlavor(h_grp);
+  H5Gclose(h_grp);
+
   /* Print the gravity parameters */
   if (e->policy & engine_policy_self_gravity) {
     h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,

src/serial_io.c

@@ -811,10 +811,18 @@ void write_output_serial(struct engine* e, const char* baseName,
                         H5P_DEFAULT);
       if (h_grp < 0) error("Error while creating SPH group");
       hydro_props_print_snapshot(h_grp, e->hydro_properties);
-      writeSPHflavour(h_grp);
+      writeChemistryFlavor(h_grp);
       H5Gclose(h_grp);
     }
 
+    /* Print the subgrid parameters */
+    h_grp = H5Gcreate(h_file, "/SubgridScheme", H5P_DEFAULT, H5P_DEFAULT,
+                      H5P_DEFAULT);
+    if (h_grp < 0) error("Error while creating subgrid group");
+    writeCoolingFlavor(h_grp);
+    writeChemistryFlavor(h_grp);
+    H5Gclose(h_grp);
+
     /* Print the gravity parameters */
     if (e->policy & engine_policy_self_gravity) {
       h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,