src/chemistry/GEAR/chemistry_struct.h

@@ -70,6 +70,10 @@ struct chemistry_bpart_data {};
 /**
  * @brief Chemical abundances traced by the #sink in the GEAR model.
  */
-struct chemistry_sink_data {};
+struct chemistry_sink_data {
+
+  /*! Total mass of element in a particle. */
+  double metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT];
+};
 
 #endif /* SWIFT_CHEMISTRY_STRUCT_GEAR_H */

src/chemistry/GEAR_DIFFUSION/chemistry.h

@@ -413,8 +413,35 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
     const struct chemistry_global_data* data, struct spart* restrict sp) {
 
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
-    sp->chemistry_data.metal_mass_fraction[i] =
-        data->initial_metallicities[i] * sp->mass;
+    /* Bug fix (26.07.2024): Check that the initial me metallicities are non
+       negative. */
+    if (data->initial_metallicities[i] >= 0) {
+      /* Use the value from the parameter file */
+      sp->chemistry_data.metal_mass_fraction[i] =
+          data->initial_metallicities[i];
+    }
+    /* else : Use the value from the IC. We are reading the metal mass
+     fraction. So do not overwrite the metallicities */
+  }
+}
+
+/**
+ * @brief Sets the chemistry properties of the sink particles to a valid start
+ * state.
+ *
+ * @param data The global chemistry information.
+ * @param sink Pointer to the sink particle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_sink(
+    const struct chemistry_global_data* data, struct sink* restrict sink) {
+
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    /* Use the value from the parameter file */
+    if (data->initial_metallicities[i] >= 0) {
+      sink->chemistry_data.metal_mass_fraction[i] =
+          data->initial_metallicities[i];
+    }
+    /* else : read the metallicities from the ICs. */
   }
 }
 

src/chemistry/GEAR_DIFFUSION/chemistry_iact.h

@@ -140,7 +140,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
 }
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (symmetric version)
  *
  * @param r2 Comoving square distance between the two particles.
@@ -226,7 +226,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_diffusion(
 }
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (nonsymmetric version)
  *
  * @param r2 Comoving square distance between the two particles.

src/chemistry/GEAR_DIFFUSION/chemistry_io.h

@@ -106,6 +106,19 @@ INLINE static int chemistry_write_sparticles(const struct spart* sparts,
   return 1;
 }
 
+/**
+ * @brief Specifies which sink fields to write to a dataset
+ *
+ * @param sinks The #sink array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_sinkparticles(const struct sink* sinks,
+                                                struct io_props* list) {
+  return 0;
+}
+
 /**
  * @brief Specifies which black hole particle fields to write to a dataset
  *

src/chemistry/QLA/chemistry.h

@@ -189,6 +189,16 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part(
 __attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
     const struct chemistry_global_data* data, struct spart* restrict sp) {}
 
+/**
+ * @brief Sets the chemistry properties of the sink particles to a valid start
+ * state.
+ *
+ * @param data The global chemistry information.
+ * @param sink Pointer to the sink particle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_sink(
+    const struct chemistry_global_data* data, struct sink* restrict sink) {}
+
 /**
  * @brief Initialise the chemistry properties of a black hole with
  * the chemistry properties of the gas it is born from.

src/chemistry/QLA/chemistry_iact.h

@@ -61,7 +61,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
     const float H) {}
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (symmetric version)
  *
  * @param r2 Comoving square distance between the two particles.
@@ -86,7 +86,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_diffusion(
     const struct cosmology *cosmo, const int with_cosmology) {}
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (nonsymmetric version)
  *
  * @param r2 Comoving square distance between the two particles.

src/chemistry/QLA/chemistry_io.h

@@ -76,6 +76,19 @@ INLINE static int chemistry_write_sparticles(const struct spart* sparts,
   return 0;
 }
 
+/**
+ * @brief Specifies which sink fields to write to a dataset
+ *
+ * @param sinks The #sink array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_sinkparticles(const struct sink* sinks,
+                                                struct io_props* list) {
+  return 0;
+}
+
 /**
  * @brief Specifies which bparticle fields to write to a dataset
  *

src/chemistry/none/chemistry.h

@@ -192,6 +192,16 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part(
 __attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
     const struct chemistry_global_data* data, struct spart* restrict sp) {}
 
+/**
+ * @brief Sets the chemistry properties of the sink particles to a valid start
+ * state.
+ *
+ * @param data The global chemistry information.
+ * @param sink Pointer to the sink particle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_sink(
+    const struct chemistry_global_data* data, struct sink* restrict sink) {}
+
 /**
  * @brief Initialise the chemistry properties of a black hole with
  * the chemistry properties of the gas it is born from.

src/chemistry/none/chemistry_iact.h

@@ -39,8 +39,9 @@
  * @param H Current Hubble parameter.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_chemistry(
-    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
-    struct part *restrict pj, const float a, const float H) {}
+    const float r2, const float dx[3], const float hi, const float hj,
+    struct part *restrict pi, struct part *restrict pj, const float a,
+    const float H) {}
 
 /**
  * @brief do chemistry computation after the runner_iact_density (non symmetric
@@ -56,11 +57,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_chemistry(
  * @param H Current Hubble parameter.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
-    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
-    const struct part *restrict pj, const float a, const float H) {}
+    const float r2, const float dx[3], const float hi, const float hj,
+    struct part *restrict pi, const struct part *restrict pj, const float a,
+    const float H) {}
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (symmetric version)
  *
  * @param r2 Comoving square distance between the two particles.
@@ -79,13 +81,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
  *
  */
 __attribute__((always_inline)) INLINE static void runner_iact_diffusion(
-    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
-    struct part *restrict pj, const float a, const float H,
-    const float time_base, const integertime_t t_current,
+    const float r2, const float dx[3], const float hi, const float hj,
+    struct part *restrict pi, struct part *restrict pj, const float a,
+    const float H, const float time_base, const integertime_t t_current,
     const struct cosmology *cosmo, const int with_cosmology) {}
 
 /**
- * @brief do metal diffusion computation in the <FORCE LOOP>
+ * @brief do metal diffusion computation in the force loop
  * (nonsymmetric version)
  *
  * @param r2 Comoving square distance between the two particles.
@@ -104,9 +106,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_diffusion(
  *
  */
 __attribute__((always_inline)) INLINE static void runner_iact_nonsym_diffusion(
-    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
-    struct part *restrict pj, const float a, const float H,
-    const float time_base, const integertime_t t_current,
+    const float r2, const float dx[3], const float hi, const float hj,
+    struct part *restrict pi, struct part *restrict pj, const float a,
+    const float H, const float time_base, const integertime_t t_current,
     const struct cosmology *cosmo, const int with_cosmology) {}
 
 #endif /* SWIFT_NONE_CHEMISTRY_IACT_H */

src/chemistry/none/chemistry_io.h

@@ -76,6 +76,19 @@ INLINE static int chemistry_write_sparticles(const struct spart* sparts,
   return 0;
 }
 
+/**
+ * @brief Specifies which sink fields to write to a dataset
+ *
+ * @param sinks The #sink array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_sinkparticles(const struct sink* sinks,
+                                                struct io_props* list) {
+  return 0;
+}
+
 /**
  * @brief Specifies which bparticle fields to write to a dataset
  *

src/common_io.c

@@ -101,6 +101,8 @@ hid_t io_hdf5_type(enum IO_DATA_TYPE type) {
       return H5T_NATIVE_DOUBLE;
     case CHAR:
       return H5T_NATIVE_CHAR;
+    case BOOL:
+      return H5T_NATIVE_HBOOL;
     default:
       error("Unknown type");
       return 0;
@@ -470,6 +472,16 @@ void io_write_attribute_i(hid_t grp, const char* name, int data) {
   io_write_attribute(grp, name, INT, &data, 1);
 }
 
+/**
+ * @brief Writes a bool value (passed as an int) as an attribute
+ * @param grp The group in which to write
+ * @param name The name of the attribute
+ * @param data The value to write
+ */
+void io_write_attribute_b(hid_t grp, const char* name, int data) {
+  io_write_attribute(grp, name, BOOL, &data, 1);
+}
+
 /**
  * @brief Writes a long value as an attribute
  * @param grp The group in which to write
@@ -1615,7 +1627,7 @@ void io_collect_gparts_neutrino_to_write(
  */
 void io_make_snapshot_subdir(const char* dirname) {
 
-  if (strcmp(dirname, ".") != 0 && strnlen(dirname, PARSER_MAX_LINE_SIZE) > 0) {
+  if (strcmp(dirname, ".") != 0 && strnlen(dirname, FILENAME_BUFFER_SIZE) > 0) {
     safe_checkdir(dirname, /*create=*/1);
   }
 }
@@ -1800,6 +1812,7 @@ void io_select_sink_fields(const struct sink* const sinks,
                            int* const num_fields, struct io_props* const list) {
 
   sink_write_particles(sinks, list, num_fields, with_cosmology);
+  *num_fields += chemistry_write_sinkparticles(sinks, list + *num_fields);
 }
 
 /**

src/common_io.h

@@ -31,6 +31,8 @@
 #define PARTICLE_GROUP_BUFFER_SIZE 50
 #define FILENAME_BUFFER_SIZE 150
 #define IO_BUFFER_ALIGNMENT 1024
+#define HDF5_LOWEST_FILE_FORMAT_VERSION H5F_LIBVER_V18
+#define HDF5_HIGHEST_FILE_FORMAT_VERSION H5F_LIBVER_LATEST
 
 /* Avoid cyclic inclusion problems */
 struct cell;
@@ -66,6 +68,7 @@ enum IO_DATA_TYPE {
   FLOAT,
   DOUBLE,
   CHAR,
+  BOOL,
   SIZE_T,
 };
 
@@ -95,6 +98,7 @@ void io_write_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
 void io_write_attribute_d(hid_t grp, const char* name, double data);
 void io_write_attribute_f(hid_t grp, const char* name, float data);
 void io_write_attribute_i(hid_t grp, const char* name, int data);
+void io_write_attribute_b(hid_t grp, const char* name, int data);
 void io_write_attribute_l(hid_t grp, const char* name, long data);
 void io_write_attribute_ll(hid_t grp, const char* name, long long data);
 void io_write_attribute_s(hid_t grp, const char* name, const char* str);

src/common_io_copy.c

@@ -488,7 +488,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
 
   } else { /* Converting particle to data */
 
-    if (props.convert_part_f != NULL) {
+    if (props.type == FLOAT && props.parts != NULL) {
 
       /* Prepare some parameters */
       float* temp_f = (float*)temp;
@@ -500,7 +500,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_part_f_mapper, temp_f, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_part_i != NULL) {
+    } else if (props.type == INT && props.parts != NULL) {
 
       /* Prepare some parameters */
       int* temp_i = (int*)temp;
@@ -512,7 +512,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_part_i_mapper, temp_i, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_part_d != NULL) {
+    } else if (props.type == DOUBLE && props.parts != NULL) {
 
       /* Prepare some parameters */
       double* temp_d = (double*)temp;
@@ -524,7 +524,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_part_d_mapper, temp_d, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_part_l != NULL) {
+    } else if (props.type == LONGLONG && props.parts != NULL) {
 
       /* Prepare some parameters */
       long long* temp_l = (long long*)temp;
@@ -536,7 +536,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_part_l_mapper, temp_l, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_gpart_f != NULL) {
+    } else if (props.type == FLOAT && props.gparts != NULL) {
 
       /* Prepare some parameters */
       float* temp_f = (float*)temp;
@@ -548,7 +548,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_gpart_f_mapper, temp_f, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_gpart_i != NULL) {
+    } else if (props.type == INT && props.gparts != NULL) {
 
       /* Prepare some parameters */
       int* temp_i = (int*)temp;
@@ -560,7 +560,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_gpart_i_mapper, temp_i, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_gpart_d != NULL) {
+    } else if (props.type == DOUBLE && props.gparts != NULL) {
 
       /* Prepare some parameters */
       double* temp_d = (double*)temp;
@@ -572,7 +572,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_gpart_d_mapper, temp_d, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_gpart_l != NULL) {
+    } else if (props.type == LONGLONG && props.gparts != NULL) {
 
       /* Prepare some parameters */
       long long* temp_l = (long long*)temp;
@@ -584,7 +584,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_gpart_l_mapper, temp_l, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_spart_f != NULL) {
+    } else if (props.type == FLOAT && props.sparts != NULL) {
 
       /* Prepare some parameters */
       float* temp_f = (float*)temp;
@@ -596,7 +596,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_spart_f_mapper, temp_f, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_spart_i != NULL) {
+    } else if (props.type == INT && props.sparts != NULL) {
 
       /* Prepare some parameters */
       int* temp_i = (int*)temp;
@@ -608,7 +608,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_spart_i_mapper, temp_i, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_spart_d != NULL) {
+    } else if (props.type == DOUBLE && props.sparts != NULL) {
 
       /* Prepare some parameters */
       double* temp_d = (double*)temp;
@@ -620,7 +620,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_spart_d_mapper, temp_d, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_spart_l != NULL) {
+    } else if (props.type == LONGLONG && props.sparts != NULL) {
 
       /* Prepare some parameters */
       long long* temp_l = (long long*)temp;
@@ -632,7 +632,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_spart_l_mapper, temp_l, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_sink_f != NULL) {
+    } else if (props.type == FLOAT && props.sinks != NULL) {
 
       /* Prepare some parameters */
       float* temp_f = (float*)temp;
@@ -644,7 +644,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_sink_f_mapper, temp_f, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_sink_i != NULL) {
+    } else if (props.type == INT && props.sinks != NULL) {
 
       /* Prepare some parameters */
       int* temp_i = (int*)temp;
@@ -656,7 +656,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_sink_i_mapper, temp_i, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_sink_d != NULL) {
+    } else if (props.type == DOUBLE && props.sinks != NULL) {
 
       /* Prepare some parameters */
       double* temp_d = (double*)temp;
@@ -668,7 +668,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_sink_d_mapper, temp_d, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_sink_l != NULL) {
+    } else if (props.type == LONGLONG && props.sinks != NULL) {
 
       /* Prepare some parameters */
       long long* temp_l = (long long*)temp;
@@ -680,7 +680,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_sink_l_mapper, temp_l, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_bpart_f != NULL) {
+    } else if (props.type == FLOAT && props.bparts != NULL) {
 
       /* Prepare some parameters */
       float* temp_f = (float*)temp;
@@ -692,7 +692,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_bpart_f_mapper, temp_f, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_bpart_i != NULL) {
+    } else if (props.type == INT && props.bparts != NULL) {
 
       /* Prepare some parameters */
       int* temp_i = (int*)temp;
@@ -704,7 +704,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_bpart_i_mapper, temp_i, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_bpart_d != NULL) {
+    } else if (props.type == DOUBLE && props.bparts != NULL) {
 
       /* Prepare some parameters */
       double* temp_d = (double*)temp;
@@ -716,7 +716,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_bpart_d_mapper, temp_d, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_bpart_l != NULL) {
+    } else if (props.type == LONGLONG && props.bparts != NULL) {
 
       /* Prepare some parameters */
       long long* temp_l = (long long*)temp;
@@ -729,7 +729,9 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      threadpool_auto_chunk_size, (void*)&props);
 
     } else {
-      error("Missing conversion function");
+
+      if (N != 0 && props.ptr_func != NULL)
+        error("Missing conversion function");
     }
   }
 

src/common_io_fields.c

@@ -400,8 +400,9 @@ void io_write_output_field_parameter(const char* filename, int with_cosmology,
         strcpy(&comment_write_buffer[PARSER_MAX_LINE_SIZE / 2 - 4], "...");
       }
 
-      fprintf(file, "  %s_%s: %s  # %s : %s\n", list[i].name,
-              part_type_names[ptype], "on", comment_write_buffer, unit_buffer);
+      fprintf(file, "  %s_%s: %s  # (%dD - %zd bytes / dim) %s : %s\n",
+              list[i].name, part_type_names[ptype], "on", list[i].dimension,
+              io_sizeof_type(list[i].type), comment_write_buffer, unit_buffer);
     }
 
     fprintf(file, "\n");

src/const.h

@@ -27,7 +27,7 @@
 
 /* Type of gradients to use (GIZMO_SPH only) */
 /* If no option is chosen, no gradients are used (first order scheme) */
-//#define GRADIENTS_SPH
+// #define GRADIENTS_SPH
 #define GRADIENTS_GIZMO
 
 /* Types of slope limiter to use (GIZMO_SPH only) */
@@ -40,11 +40,11 @@
 
 /* Options to control the movement of particles for GIZMO_SPH. */
 /* This option disables particle movement */
-//#define GIZMO_FIX_PARTICLES
+// #define GIZMO_FIX_PARTICLES
 /* Try to keep cells regular by adding a correction velocity. */
-//#define GIZMO_STEER_MOTION
+// #define GIZMO_STEER_MOTION
 /* Use the total energy instead of the thermal energy as conserved variable. */
-//#define GIZMO_TOTAL_ENERGY
+// #define GIZMO_TOTAL_ENERGY
 
 /* Options to control handling of unphysical values (GIZMO_SPH only). */
 /* In GIZMO, mass and energy (and hence density and pressure) can in principle
@@ -53,7 +53,7 @@
    If no options are selected below, we assume (and pray) that this will not
    happen, and add no restrictions to how these variables are treated. */
 /* Check for unphysical values and crash if they occur. */
-//#define GIZMO_UNPHYSICAL_ERROR
+// #define GIZMO_UNPHYSICAL_ERROR
 /* Check for unphysical values and reset them to safe values. */
 #define GIZMO_UNPHYSICAL_RESCUE
 /* Show a warning message if an unphysical value was reset (only works if
@@ -65,9 +65,9 @@
 /* Parameters that control how GIZMO handles pathological particle
    configurations. */
 /* Show a warning message if a pathological configuration has been detected. */
-//#define GIZMO_PATHOLOGICAL_WARNING
+// #define GIZMO_PATHOLOGICAL_WARNING
 /* Crash if a pathological configuration has been detected. */
-//#define GIZMO_PATHOLOGICAL_ERROR
+// #define GIZMO_PATHOLOGICAL_ERROR
 /* Maximum allowed gradient matrix condition number. If the condition number of
    the gradient matrix (defined in equation C1 in Hopkins, 2015) is larger than
    this value, we artificially increase the number of neighbours to get a more
@@ -81,26 +81,14 @@
    gradient matrix and use SPH gradients instead. */
 #define const_gizmo_min_wcorr 0.5f
 
-/* Types of gradients to use for SHADOWFAX_SPH */
-/* If no option is chosen, no gradients are used (first order scheme) */
-#define SHADOWFAX_GRADIENTS
-
-/* SHADOWFAX_SPH slope limiters */
-#define SHADOWFAX_SLOPE_LIMITER_PER_FACE
-#define SHADOWFAX_SLOPE_LIMITER_CELL_WIDE
-
-/* Options to control SHADOWFAX_SPH */
-/* This option disables cell movement */
-//#define SHADOWFAX_FIX_CELLS
-/* This option enables cell steering, i.e. trying to keep the cells regular by
-   adding a correction to the cell velocities.*/
-#define SHADOWFAX_STEER_CELL_MOTION
-/* This option evolves the total energy instead of the thermal energy */
-//#define SHADOWFAX_TOTAL_ENERGY
+/* Options controlling ShadowSWIFT */
+/* Options controlling acceleration strategies*/
+/*! @brief Option enabling a more relaxed completeness criterion */
+#define SHADOWSWIFT_RELAXED_COMPLETENESS
 
 /* Source terms */
 #define SOURCETERMS_NONE
-//#define SOURCETERMS_SN_FEEDBACK
+// #define SOURCETERMS_SN_FEEDBACK
 
 /* GRACKLE doesn't really like exact zeroes, so use something
  * comparatively small instead. */

src/cooling/EAGLE/cooling.c

@@ -693,6 +693,30 @@ float cooling_get_temperature(
   return exp10(log_10_T);
 }
 
+/**
+ * @brief Compute the electron number density of a #part based on the cooling
+ * function.
+ *
+ * Does not exist in this model. We return 0.
+ *
+ * @param phys_const #phys_const data structure.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param us The internal system of units.
+ * @param cosmo #cosmology data structure.
+ * @param cooling #cooling_function_data struct.
+ * @param p #part data.
+ * @param xp Pointer to the #xpart data.
+ */
+double cooling_get_electron_density(const struct phys_const *phys_const,
+                                    const struct hydro_props *hydro_props,
+                                    const struct unit_system *us,
+                                    const struct cosmology *cosmo,
+                                    const struct cooling_function_data *cooling,
+                                    const struct part *p,
+                                    const struct xpart *xp) {
+  return 0.;
+}
+
 /**
  * @brief Compute the electron pressure of a #part based on the cooling
  * function.

src/cooling/EAGLE/cooling.h

@@ -60,6 +60,14 @@ float cooling_timestep(const struct cooling_function_data *cooling,
                        const struct hydro_props *hydro_props,
                        const struct part *p, const struct xpart *xp);
 
+double cooling_get_electron_density(const struct phys_const *phys_const,
+                                    const struct hydro_props *hydro_props,
+                                    const struct unit_system *us,
+                                    const struct cosmology *cosmo,
+                                    const struct cooling_function_data *cooling,
+                                    const struct part *p,
+                                    const struct xpart *xp);
+
 double cooling_get_electron_pressure(
     const struct phys_const *phys_const, const struct hydro_props *hydro_props,
     const struct unit_system *us, const struct cosmology *cosmo,
@@ -170,7 +178,4 @@ void cooling_print_backend(const struct cooling_function_data *cooling);
 
 void cooling_clean(struct cooling_function_data *data);
 
-/*! Stub defined to let the BH model compile */
-#define colibre_cooling_N_elementtypes 1
-
 #endif /* SWIFT_COOLING_EAGLE_H */

src/cooling/PS2020/cooling_io.h

@@ -177,18 +177,18 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
       "Temperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts, xparts,
       convert_part_T, "Temperatures of the gas particles");
 
-  list[1] = io_make_output_field_convert_part(
+  list[1] = io_make_physical_output_field_convert_part(
       "SubgridTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts,
-      xparts, convert_part_sub_T,
+      xparts, /*can convert to comoving=*/1, convert_part_sub_T,
       "The subgrid temperatures if the particles are within deltaT of the "
       "entropy floor the subgrid temperature is calculated assuming a "
       "pressure equilibrium on the entropy floor, if the particles are "
       "above deltaT of the entropy floor the subgrid temperature is "
       "identical to the SPH temperature.");
 
-  list[2] = io_make_output_field_convert_part(
-      "SubgridPhysicalDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, parts,
-      xparts, convert_part_sub_rho,
+  list[2] = io_make_physical_output_field_convert_part(
+      "SubgridPhysicalDensities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f, parts,
+      xparts, /*can convert to comoving=*/1, convert_part_sub_rho,
       "The subgrid physical density if the particles are within deltaT of the "
       "entropy floor the subgrid density is calculated assuming a pressure "
       "equilibrium on the entropy floor, if the particles are above deltaT "
@@ -222,15 +222,15 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
       "floor, by extrapolating to the equilibrium curve assuming constant "
       "pressure.");
 
-  list[6] = io_make_output_field_convert_part(
-      "ElectronNumberDensities", DOUBLE, 1, UNIT_CONV_NUMBER_DENSITY, 0.f,
-      parts, xparts, convert_part_e_density,
+  list[6] = io_make_physical_output_field_convert_part(
+      "ElectronNumberDensities", DOUBLE, 1, UNIT_CONV_NUMBER_DENSITY, -3.f,
+      parts, xparts, /*can convert to comoving=*/1, convert_part_e_density,
       "Electron number densities in the physical frame computed based on the "
       "cooling tables. This is 0 for star-forming particles.");
 
-  list[7] = io_make_output_field_convert_part(
+  list[7] = io_make_physical_output_field_convert_part(
       "ComptonYParameters", DOUBLE, 1, UNIT_CONV_AREA, 0.f, parts, xparts,
-      convert_part_y_compton,
+      /*can convert to comoving=*/0, convert_part_y_compton,
       "Compton y parameters in the physical frame computed based on the "
       "cooling tables. This is 0 for star-forming particles.");
 

src/cooling/const_lambda/cooling_io.h

@@ -80,9 +80,9 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
       "Temperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts, xparts,
       convert_part_T, "Temperatures of the gas particles");
 
-  list[1] = io_make_output_field(
+  list[1] = io_make_physical_output_field(
       "RadiatedEnergies", FLOAT, 1, UNIT_CONV_ENERGY, 0.f, xparts,
-      cooling_data.radiated_energy,
+      cooling_data.radiated_energy, /*can convert to comoving=*/0,
       "Thermal energies radiated by the cooling mechanism");
 
   return 2;

src/cooling/grackle/cooling.c

@@ -64,6 +64,22 @@ gr_float cooling_time(const struct phys_const* phys_const,
                       const struct cooling_function_data* cooling,
                       const struct part* p, struct xpart* xp);
 
+double cooling_get_physical_density(
+    const struct part* p, const struct cosmology* cosmo,
+    const struct cooling_function_data* cooling);
+
+/**
+ * @brief Record the time when cooling was switched off for a particle.
+ *
+ * @param p #part data.
+ * @param xp Pointer to the #xpart data.
+ * @param time The time when the cooling was switched off.
+ */
+INLINE void cooling_set_part_time_cooling_off(struct part* p, struct xpart* xp,
+                                              const double time) {
+  xp->cooling_data.time_last_event = time;
+}
+
 /**
  * @brief Common operations performed on the cooling function at a
  * given time-step or redshift.
@@ -207,32 +223,166 @@ void cooling_first_init_part(const struct phys_const* phys_const,
 #if COOLING_GRACKLE_MODE >= 1
   gr_float zero = 1.e-20;
 
-  /* NOTE: at this stage, we assume neutral gas
-   * a better determination will be done in cooling_post_init_part */
+  /* NOTE: if the ratio with respect to hydrogen is given, use it. Instead,
+   * we assume neutral gas.
+   * A better determination of the abundances can be done in
+   * cooling_post_init_part */
+
+  /* Compute nH (formally divided by the gas density and assuming the proton
+   * mass to be one) */
+  double nH = 1.0;
+  if (grackle_data != NULL) {
+    nH = grackle_data->HydrogenFractionByMass;
+  }
+
+  /* Compute nHe (formally divided by the gas density and assuming the proton
+   * mass to be one) */
+  double nHe = 0.0;
+  if (grackle_data != NULL) {
+    nHe = (1.f - grackle_data->HydrogenFractionByMass) / 4.f;
+  }
+
+  /* Electron density */
+  double ne = zero;
+
+  /* Sum of mass fraction (to test its consistency) */
+  double Xtot = 0.f;
 
   /* primordial chemistry >= 1 */
-  /* assume neutral gas */
-  xp->cooling_data.HI_frac = grackle_data->HydrogenFractionByMass;
-  xp->cooling_data.HII_frac = zero;
-  xp->cooling_data.HeI_frac = 1. - grackle_data->HydrogenFractionByMass;
-  xp->cooling_data.HeII_frac = zero;
-  xp->cooling_data.HeIII_frac = zero;
-  xp->cooling_data.e_frac = zero;
+
+  /* Hydrogen I */
+  if (cooling->initial_nHII_to_nH_ratio >= 0.f)
+    xp->cooling_data.HI_frac = nH * (1.f - cooling->initial_nHII_to_nH_ratio);
+  else
+    xp->cooling_data.HI_frac = nH;
+
+  Xtot += xp->cooling_data.HI_frac;
+
+  /* Hydrogen II */
+  if (cooling->initial_nHII_to_nH_ratio >= 0.f) {
+    double nHII = nH * cooling->initial_nHII_to_nH_ratio;
+    xp->cooling_data.HII_frac = nHII;
+    ne += nHII;
+  } else
+    xp->cooling_data.HII_frac = zero;
+
+  Xtot += xp->cooling_data.HII_frac;
+
+  /* Helium I */
+  if (cooling->initial_nHeI_to_nH_ratio >= 0.f) {
+    double nHeI = nH * cooling->initial_nHeI_to_nH_ratio;
+    xp->cooling_data.HeI_frac = nHeI * 4.f;
+  } else
+    xp->cooling_data.HeI_frac = nHe * 4.f;
+
+  Xtot += xp->cooling_data.HeI_frac;
+
+  /* Helium II */
+  if (cooling->initial_nHeII_to_nH_ratio >= 0.f) {
+    double nHeII = nH * cooling->initial_nHeII_to_nH_ratio;
+    xp->cooling_data.HeII_frac = nHeII * 4.f;
+    ne += nHeII;
+  } else
+    xp->cooling_data.HeII_frac = zero;
+
+  Xtot += xp->cooling_data.HeII_frac;
+
+  /* Helium III */
+  if (cooling->initial_nHeIII_to_nH_ratio >= 0.f) {
+    double nHeIII = nH * cooling->initial_nHeIII_to_nH_ratio;
+    xp->cooling_data.HeIII_frac = nHeIII * 4.f;
+    ne += 2.f * nHeIII;
+  } else
+    xp->cooling_data.HeIII_frac = zero;
+
+  Xtot += xp->cooling_data.HeIII_frac;
+
+  /* electron mass fraction (multiplied by the proton mass (Grackle convention)
+   */
+  xp->cooling_data.e_frac = ne;
+
 #endif  // MODE >= 1
 
 #if COOLING_GRACKLE_MODE >= 2
   /* primordial chemistry >= 2 */
-  xp->cooling_data.HM_frac = zero;
-  xp->cooling_data.H2I_frac = zero;
-  xp->cooling_data.H2II_frac = zero;
+
+  /* Hydrogen- */
+  if (cooling->initial_nHM_to_nH_ratio >= 0.f) {
+    double nHM = nH * cooling->initial_nHM_to_nH_ratio;
+    xp->cooling_data.HM_frac = nHM;
+    ne -= nHM;
+  } else
+    xp->cooling_data.HM_frac = zero;
+
+  Xtot += xp->cooling_data.HM_frac;
+
+  /* H2I */
+  if (cooling->initial_nH2I_to_nH_ratio >= 0.f) {
+    double nH2I = nH * cooling->initial_nH2I_to_nH_ratio;
+    xp->cooling_data.H2I_frac = nH2I * 2.f;
+  } else
+    xp->cooling_data.H2I_frac = zero;
+
+  Xtot += xp->cooling_data.H2I_frac;
+
+  /* H2II */
+  if (cooling->initial_nH2II_to_nH_ratio >= 0.f) {
+    double nH2II = nH * cooling->initial_nH2II_to_nH_ratio;
+    xp->cooling_data.H2II_frac = nH2II * 2.f;
+    ne += nH2II;
+  } else
+    xp->cooling_data.H2II_frac = zero;
+
+  Xtot += xp->cooling_data.H2II_frac;
+
+  /* electron mass fraction (multiplied by the proton mass (Grackle convention)
+   */
+  xp->cooling_data.e_frac = ne;
+
 #endif  // MODE >= 2
 
 #if COOLING_GRACKLE_MODE >= 3
   /* primordial chemistry >= 3 */
-  xp->cooling_data.DI_frac = grackle_data->DeuteriumToHydrogenRatio *
-                             grackle_data->HydrogenFractionByMass;
-  xp->cooling_data.DII_frac = zero;
-  xp->cooling_data.HDI_frac = zero;
+
+  /* Deuterium I */
+  if (cooling->initial_nDI_to_nH_ratio >= 0.f) {
+    double nDI = nH * cooling->initial_nDI_to_nH_ratio;
+    xp->cooling_data.DI_frac = nDI * 2.f;
+  } else {
+    if (grackle_data != NULL) {
+      xp->cooling_data.DI_frac = grackle_data->DeuteriumToHydrogenRatio *
+                                 grackle_data->HydrogenFractionByMass;
+    }
+  }
+
+  Xtot += xp->cooling_data.DI_frac;
+
+  /* Deuterium II */
+  if (cooling->initial_nDII_to_nH_ratio >= 0.f) {
+    double nDII = nH * cooling->initial_nDII_to_nH_ratio;
+    xp->cooling_data.DII_frac = nDII * 2.f;
+    ne += nDII;
+  } else
+    xp->cooling_data.DII_frac = zero;
+
+  Xtot += xp->cooling_data.DII_frac;
+
+  /* HD I */
+  if (cooling->initial_nHDI_to_nH_ratio >= 0.f) {
+    double nHDI = nH * cooling->initial_nHDI_to_nH_ratio;
+    xp->cooling_data.HDI_frac = nHDI * 3.f;
+  } else
+    xp->cooling_data.HDI_frac = zero;
+
+  Xtot += xp->cooling_data.HDI_frac;
+
+  /* electron mass fraction (multiplied by the proton mass (Grackle convention)
+   */
+  xp->cooling_data.e_frac = ne;
+
+  if (fabs(Xtot - 1.0) > 1e-3)
+    error("Got total mass fraction of gas = %.6g", Xtot);
+
 #endif  // MODE >= 3
 }
 
@@ -284,6 +434,10 @@ float cooling_get_radiated_energy(const struct xpart* xp) {
  */
 void cooling_print_backend(const struct cooling_function_data* cooling) {
 
+  if (engine_rank != 0) {
+    return;
+  }
+
   message("Cooling function is 'Grackle'.");
   message("Using Grackle = %i", cooling->chemistry_data.use_grackle);
   message("Chemical network = %i",
@@ -293,6 +447,7 @@ void cooling_print_backend(const struct cooling_function_data* cooling) {
   message("UV background = %d", cooling->with_uv_background);
   message("Metal cooling = %i", cooling->chemistry_data.metal_cooling);
   message("Self Shielding = %i", cooling->self_shielding_method);
+  message("Maximal density = %e", cooling->cooling_density_max);
   if (cooling->self_shielding_method == -1) {
     message("Self Shelding density = %g", cooling->self_shielding_threshold);
   }
@@ -311,7 +466,26 @@ void cooling_print_backend(const struct cooling_function_data* cooling) {
           cooling->RT_HeII_ionization_rate);
   message("grackle_chemistry_data.RT_H2_dissociation_rate = %g",
           cooling->RT_H2_dissociation_rate);
-
+  message("cooling.initial_nHII_to_nH_ratio= %g",
+          cooling->initial_nHII_to_nH_ratio);
+  message("cooling.initial_nHeI_to_nH_ratio= %g",
+          cooling->initial_nHeI_to_nH_ratio);
+  message("cooling.initial_nHeII_to_nH_ratio= %g",
+          cooling->initial_nHeII_to_nH_ratio);
+  message("cooling.initial_nHeIII_to_nH_ratio= %g",
+          cooling->initial_nHeIII_to_nH_ratio);
+  message("cooling.initial_nDI_to_nH_ratio= %g",
+          cooling->initial_nDI_to_nH_ratio);
+  message("cooling.initial_nDII_to_nH_ratio= %g",
+          cooling->initial_nDII_to_nH_ratio);
+  message("cooling.initial_nHM_to_nH_ratio= %g",
+          cooling->initial_nHM_to_nH_ratio);
+  message("cooling.initial_nH2I_to_nH_ratio= %g",
+          cooling->initial_nH2I_to_nH_ratio);
+  message("cooling.initial_nH2II_to_nH_ratio= %g",
+          cooling->initial_nH2II_to_nH_ratio);
+  message("cooling.initial_nHDI_to_nH_ratio= %g",
+          cooling->initial_nHDI_to_nH_ratio);
   message("Units:");
   message("\tComoving = %i", cooling->units.comoving_coordinates);
   message("\tLength = %g", cooling->units.length_units);
@@ -770,7 +944,7 @@ gr_float cooling_new_energy(const struct phys_const* phys_const,
   data.grid_end = grid_end;
 
   /* general particle data */
-  gr_float density = hydro_get_physical_density(p, cosmo);
+  gr_float density = cooling_get_physical_density(p, cosmo, cooling);
   gr_float energy = hydro_get_physical_internal_energy(p, xp, cosmo) +
                     dt_therm * hydro_get_physical_internal_energy_dt(p, cosmo);
   energy = max(energy, hydro_props->minimal_internal_energy);
@@ -848,7 +1022,7 @@ gr_float cooling_time(const struct phys_const* phys_const,
   data.grid_end = grid_end;
 
   /* general particle data */
-  gr_float density = hydro_get_physical_density(p, cosmo);
+  gr_float density = cooling_get_physical_density(p, cosmo, cooling);
   gr_float energy = hydro_get_physical_internal_energy(p, xp, cosmo);
   energy = max(energy, hydro_props->minimal_internal_energy);
 
@@ -1106,8 +1280,10 @@ void cooling_init_units(const struct unit_system* us,
 void cooling_init_grackle(struct cooling_function_data* cooling) {
 
 #ifdef SWIFT_DEBUG_CHECKS
-  /* enable verbose for grackle */
-  grackle_verbose = 1;
+  /* Enable verbose for grackle for rank 0 only. */
+  if (engine_rank == 0) {
+    grackle_verbose = 1;
+  }
 #endif
 
   chemistry_data* chemistry = &cooling->chemistry_data;
@@ -1195,7 +1371,7 @@ void cooling_init_backend(struct swift_params* parameter_file,
     error("Grackle with multiple particles not implemented");
 
   /* read parameters */
-  cooling_read_parameters(parameter_file, cooling, phys_const);
+  cooling_read_parameters(parameter_file, cooling, phys_const, us);
 
   /* Set up the units system. */
   cooling_init_units(us, phys_const, cooling);
@@ -1247,3 +1423,41 @@ void cooling_struct_restore(struct cooling_function_data* cooling, FILE* stream,
   /* Set up grackle */
   cooling_init_grackle(cooling);
 }
+
+/**
+ * @brief Get the density of the #part. If the density is bigger than
+ * cooling_density_max, then we floor the density to this value.
+ *
+ * This function ensures that we pass to grackle a density value that is not
+ * to big to ensure good working of grackle.
+ *
+ * Note: This function is called in cooling_time() and cooling_new_energy().
+ *
+ * @param p #part data.
+ * @param cosmo #cosmology data structure.
+ * @param cooling #cooling_function_data struct.
+ */
+double cooling_get_physical_density(
+    const struct part* p, const struct cosmology* cosmo,
+    const struct cooling_function_data* cooling) {
+
+  const double part_density = hydro_get_physical_density(p, cosmo);
+  const double cooling_max_density = cooling->cooling_density_max;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (cooling_max_density > 0.0 && part_density > cooling_max_density) {
+    warning(
+        "Gas particle %lld physical density (%e) is higher than the maximal "
+        "physical density set in the parameter files (%e).",
+        p->id, part_density, cooling_max_density);
+  }
+#endif
+  /* Maximal density cooling is defined */
+  if (cooling_max_density > 0.0) {
+    return fminf(part_density, cooling_max_density);
+  }
+  /* else ( if cooling_max_density <= 0) we do not want to use a density
+     threshold, then return the part density. */
+
+  return part_density;
+}