bugfix and cleaned code

configure.ac

@@ -533,7 +533,7 @@ if test "x$with_grackle" != "xno"; then
       [grackle],
       [initialize_chemistry_data],
       [AC_DEFINE([HAVE_GRACKLE],1,[The GRACKLE library appears to be present.])
-        AC_DEFINE([CONFIG_BFLOAT_4],1,[Use double in grackle])
+        AC_DEFINE([CONFIG_BFLOAT_8],1,[Use double in grackle])
       ],
       [AC_MSG_ERROR(Cannot find grackle library!)],
       [$GRACKLE_LIBS $GRACKLE_INCS $FCLIBS]

src/cooling/grackle/cooling.h

@@ -267,6 +267,180 @@ __attribute__((always_inline)) INLINE static void cooling_print_backend(
 #endif
 }
 
+/**
+ * @brief allocate required field
+ *
+ * @param data the #grackle_field_data
+ */
+__attribute__((always_inline)) INLINE static void cooling_malloc_data(grackle_field_data *data) {
+
+#if COOLING_GRACKLE_MODE >= 1
+  /* primordial chemistry >= 1 */
+  data->HI_density = malloc(sizeof(gr_float));
+  data->HII_density = malloc(sizeof(gr_float));
+  data->HeI_density = malloc(sizeof(gr_float));
+  data->HeII_density = malloc(sizeof(gr_float));
+  data->HeIII_density = malloc(sizeof(gr_float));
+  data->e_density = malloc(sizeof(gr_float));
+#endif // MODE >= 1
+
+#if COOLING_GRACKLE_MODE >= 2
+  /* primordial chemistry >= 2 */
+  data->HM_density = malloc(sizeof(gr_float));
+  data->H2I_density = malloc(sizeof(gr_float));
+  data->H2II_density = malloc(sizeof(gr_float));
+#endif // MODE 2
+
+#if COOLING_GRACKLE_MODE >= 3
+  /* primordial chemistry >= 3 */
+  data->DI_density = malloc(sizeof(gr_float));
+  data->DII_density = malloc(sizeof(gr_float));
+  data->HDI_density = malloc(sizeof(gr_float));
+#endif // MODE >= 3
+
+  /* metal cooling = 1 */
+  data->metal_density = malloc(sizeof(gr_float));
+
+  /* /\* volumetric heating rate *\/ */
+  /* data->volumetric_heating_rate = NULL; */
+
+  /* /\* specific heating rate *\/ */
+  /* data->specific_heating_rate = NULL; */
+
+}
+
+/**
+ * @brief free the allocated memory
+ *
+ * @param data the #grackle_field_data
+ */
+
+__attribute__((always_inline)) INLINE static void cooling_free_data(grackle_field_data *data) {
+
+#if COOLING_GRACKLE_MODE >= 1
+  /* primordial chemistry >= 1 */
+  free(data->HI_density);
+  free(data->HII_density);
+  free(data->HeI_density);
+  free(data->HeII_density);
+  free(data->HeIII_density);
+  free(data->e_density);
+#endif // MODE >= 1
+
+#if COOLING_GRACKLE_MODE >= 2
+  /* primordial chemistry >= 2 */
+  free(data->HM_density);
+  free(data->H2I_density);
+  free(data->H2II_density);
+#endif // MODE 2
+
+#if COOLING_GRACKLE_MODE >= 3
+  /* primordial chemistry >= 3 */
+  free(data->DI_density);
+  free(data->DII_density);
+  free(data->HDI_density);
+#endif // MODE >= 3
+
+  /* metal cooling = 1 */
+  free(data->metal_density);
+
+  /* /\* volumetric heating rate *\/ */
+  /* data->volumetric_heating_rate = NULL; */
+
+  /* /\* specific heating rate *\/ */
+  /* data->specific_heating_rate = NULL; */
+
+}
+
+/**
+ * @brief copy xp to data
+ *
+ * requires the particle to have been transformed into density
+ *
+ * @param data the #grackle_field_data
+ * @param xp the #xpart
+ */
+__attribute__((always_inline)) INLINE static void cooling_copy_to_data(grackle_field_data *data, const struct xpart *xp) {
+
+#if COOLING_GRACKLE_MODE >= 1
+  /* primordial chemistry >= 1 */
+  data->HI_density[0] = xp->cooling_data.HI_frac;
+  data->HII_density[0] = xp->cooling_data.HII_frac;
+  data->HeI_density[0] = xp->cooling_data.HeI_frac;
+  data->HeII_density[0] = xp->cooling_data.HeII_frac;
+  data->HeIII_density[0] = xp->cooling_data.HeIII_frac;
+  data->e_density[0] = xp->cooling_data.e_frac;
+#endif // MODE >= 1
+
+#if COOLING_GRACKLE_MODE >= 2
+  /* primordial chemistry >= 2 */
+  data->HM_density[0] = xp->cooling_data.HM_frac;
+  data->H2I_density[0] = xp->cooling_data.H2I_frac;
+  data->H2II_density[0] = xp->cooling_data.H2II_frac;
+#endif // MODE 2
+
+#if COOLING_GRACKLE_MODE >= 3
+  /* primordial chemistry >= 3 */
+  data->DI_density[0] = xp->cooling_data.DI_frac;
+  data->DII_density[0] = xp->cooling_data.DII_frac;
+  data->HDI_density[0] = xp->cooling_data.HDI_frac;
+#endif // MODE >= 3
+
+  /* metal cooling = 1 */
+  data->metal_density[0] = xp->cooling_data.metal_frac;
+
+  /* volumetric heating rate */
+  data->volumetric_heating_rate = NULL;
+
+  /* specific heating rate */
+  data->specific_heating_rate = NULL;
+}
+
+
+/**
+ * @brief copy data to xp
+ *
+ * @param data the #grackle_field_data
+ * @param xp the #xpart
+ */
+__attribute__((always_inline)) INLINE static void cooling_copy_to_particle(const grackle_field_data *data, struct xpart *xp) {
+
+#if COOLING_GRACKLE_MODE >= 1
+  /* primordial chemistry >= 1 */
+  xp->cooling_data.HI_frac = data->HI_density[0];
+  xp->cooling_data.HII_frac = data->HII_density[0];
+  xp->cooling_data.HeI_frac = data->HeI_density[0];
+  xp->cooling_data.HeII_frac = data->HeII_density[0];
+  xp->cooling_data.HeIII_frac = data->HeIII_density[0];
+  xp->cooling_data.e_frac = data->e_density[0];
+#endif // MODE >= 1
+
+#if COOLING_GRACKLE_MODE >= 2
+  /* primordial chemistry >= 2 */
+  xp->cooling_data.HM_frac = data->HM_density[0];
+  xp->cooling_data.H2I_frac = data->H2I_density[0];
+  xp->cooling_data.H2II_frac = data->H2II_density[0];
+#endif // MODE 2
+
+#if COOLING_GRACKLE_MODE >= 3
+  /* primordial chemistry >= 3 */
+  xp->cooling_data.DI_frac = data->DI_density[0];
+  xp->cooling_data.DII_frac = data->DII_density[0];
+  xp->cooling_data.HDI_frac = data->HDI_density[0];
+#endif // MODE >= 3
+
+  /* metal cooling = 1 */
+  xp->cooling_data.metal_frac = data->metal_density[0];
+
+  /* /\* volumetric heating rate *\/ */
+  /* data->volumetric_heating_rate = NULL; */
+
+  /* /\* specific heating rate *\/ */
+  /* data->specific_heating_rate = NULL; */
+
+}
+
+
 /**
  * @brief Compute the cooling rate and update the particle chemistry data
  *
@@ -327,48 +501,26 @@ __attribute__((always_inline)) INLINE static double cooling_rate(
   /* transform gas fraction to densities */
   cooling_compute_density(xp, *data.density);
 
-#if COOLING_GRACKLE_MODE >= 1
-  /* primordial chemistry >= 1 */
-  data.HI_density = &xp->cooling_data.HI_frac;
-  data.HII_density = &xp->cooling_data.HII_frac;
-  data.HeI_density = &xp->cooling_data.HeI_frac;
-  data.HeII_density = &xp->cooling_data.HeII_frac;
-  data.HeIII_density = &xp->cooling_data.HeIII_frac;
-  data.e_density = &xp->cooling_data.e_frac;
-#endif // MODE >= 1
-
-#if COOLING_GRACKLE_MODE >= 2
-  /* primordial chemistry >= 2 */
-  data.HM_density = &xp->cooling_data.HM_frac;
-  data.H2I_density = &xp->cooling_data.H2I_frac;
-  data.H2II_density = &xp->cooling_data.H2II_frac;
-#endif // MODE 2
-
-#if COOLING_GRACKLE_MODE >= 3
-  /* primordial chemistry >= 3 */
-  data.DI_density = &xp->cooling_data.DI_frac;
-  data.DII_density = &xp->cooling_data.DII_frac;
-  data.HDI_density = &xp->cooling_data.HDI_frac;
-#endif // MODE >= 3
-
+  /* allocate grackle data */
+  cooling_malloc_data(&data);
 
-  /* metal cooling = 1 */
-  data.metal_density = &xp->cooling_data.metal_frac;
-
-  /* volumetric heating rate */
-  data.volumetric_heating_rate = NULL;
-
-  /* specific heating rate */
-  data.specific_heating_rate = NULL;
+  /* copy data from particle to grackle data */
+  cooling_copy_to_data(&data, xp);
 
   /* solve chemistry with table */
   if (solve_chemistry(&units, &data, dt) == 0) {
     error("Error in solve_chemistry.");
   }
 
+  /* copy from grackle data to particle */
+  cooling_copy_to_particle(&data, xp);
+
   /* transform densities to gas fraction */
   cooling_compute_fraction(xp, *data.density);
 
+  /* free allocated memory */
+  cooling_free_data(&data);
+
   /* compute rate */
   return (energy - energy_before) / dt;
 }

src/cooling/grackle/cooling_struct.h

@@ -61,24 +61,24 @@ struct cooling_xpart_data {
   /* here all fractions are mass fraction */
 #if COOLING_GRACKLE_MODE >= 1
   /* primordial chemistry >= 1 */
-  gr_float HI_frac;
-  gr_float HII_frac;
-  gr_float HeI_frac;
-  gr_float HeII_frac;
-  gr_float HeIII_frac;
-  gr_float e_frac;
+  float HI_frac;
+  float HII_frac;
+  float HeI_frac;
+  float HeII_frac;
+  float HeIII_frac;
+  float e_frac;
 
 #if COOLING_GRACKLE_MODE >= 2
   /* primordial chemistry >= 2 */
-  gr_float HM_frac;
-  gr_float H2I_frac;
-  gr_float H2II_frac;
+  float HM_frac;
+  float H2I_frac;
+  float H2II_frac;
 
 #if COOLING_GRACKLE_MODE >= 3
   /* primordial chemistry >= 3 */
-  gr_float DI_frac;
-  gr_float DII_frac;
-  gr_float HDI_frac;
+  float DI_frac;
+  float DII_frac;
+  float HDI_frac;
 #endif // MODE >= 3
 
 #endif // MODE >= 2
@@ -86,7 +86,7 @@ struct cooling_xpart_data {
 #endif // MODE >= 1
   
   /* metal cooling = 1 */
-  gr_float metal_frac;
+  float metal_frac;
 };
 
 #endif /* SWIFT_COOLING_STRUCT_NONE_H */