Moved init from wrapper to cooling, changed a few variables and cleaned the code

src/cooling/grackle/cooling.h

@@ -27,6 +27,7 @@
 /* Some standard headers. */
 #include <float.h>
 #include <math.h>
+#include <grackle.h>
 
 /* Local includes. */
 #include "error.h"
@@ -36,13 +37,37 @@
 #include "physical_constants.h"
 #include "units.h"
 
-/* include the grackle wrapper */
-#include "grackle_wrapper.h"
+/* need to rework code if changed */
+#define GRACKLE_NPART 1
+
 
 /**
- * @brief Compute the cooling rate
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state.
  *
- * We do nothing.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void cooling_init_part(
+    const struct part* restrict p, struct xpart* restrict xp) {
+
+  xp->cooling_data.radiated_energy = 0.f;
+}
+
+/**
+ * @brief Returns the total radiated energy by this particle.
+ *
+ * @param xp The extended particle data
+ */
+__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
+    const struct xpart* restrict xp) {
+
+  return xp->cooling_data.radiated_energy;
+}
+
+
+/**
+ * @brief Compute the cooling rate
  *
  * @param phys_const The physical constants in internal units.
  * @param us The internal system of units.
@@ -134,30 +159,6 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
   return FLT_MAX;
 }
 
-/**
- * @brief Sets the cooling properties of the (x-)particles to a valid start
- * state.
- *
- * @param p Pointer to the particle data.
- * @param xp Pointer to the extended particle data.
- */
-__attribute__((always_inline)) INLINE static void cooling_init_part(
-    const struct part* restrict p, struct xpart* restrict xp) {
-
-  xp->cooling_data.radiated_energy = 0.f;
-}
-
-/**
- * @brief Returns the total radiated energy by this particle.
- *
- * @param xp The extended particle data
- */
-__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
-    const struct xpart* restrict xp) {
-
-  return xp->cooling_data.radiated_energy;
-}
-
 /**
  * @brief Initialises the cooling properties.
  *
@@ -171,27 +172,67 @@ static INLINE void cooling_init_backend(
     const struct phys_const* phys_const,
     struct cooling_function_data* cooling) {
 
-  double units_density, units_length, units_time;
-  int grackle_chemistry;
-  int UVbackground;
-
+    /* read parameters */
   parser_get_param_string(parameter_file, "GrackleCooling:GrackleCloudyTable",
-                          cooling->GrackleCloudyTable);
-  cooling->UVbackground =
-      parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");
-  cooling->GrackleRedshift =
-      parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");
-  cooling->GrackleHSShieldingDensityThreshold = parser_get_param_double(
+                          cooling->cloudy_table);
+  cooling->uv_background =
+    parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");
+
+  cooling->redshift =
+    parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");
+
+  cooling->density_self_shielding = parser_get_param_double(
       parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold");
+  
+#ifdef SWIFT_DEBUG_CHECKS
+  /* enable verbose for grackle */
+  grackle_verbose = 1;
+#endif
 
-  UVbackground = cooling->UVbackground;
-  grackle_chemistry = 0; /* forced to be zero : read table */
+  /* Set up the units system.
+     These are conversions from code units to cgs. */
+
+  /* first cosmo */
+  cooling->units.a_units = 1.0;  // units for the expansion factor (1/1+zi)
+
+  /* We assume here all physical quantities to
+     be in proper coordinate (not comobile)  */
+  cooling->comoving_coordinates = 0;
+
+  /* then units */
+  cooling->units.density_units = us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
+  cooling->units.length_units = us->UnitLength_in_cgs;
+  cooling->units.time_units = us->UnitTime_in_cgs;
+  cooling->units.velocity_units =
+    cooling->units.a_units * cooling->units.length_units / cooling->units.time_units;
+ 
+  /* Create a chemistry object for parameters and rate data. */
+  if (set_default_chemistry_parameters() == 0) {
+    error("Error in set_default_chemistry_parameters.");
+  }
 
-  units_density = us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
-  units_length = us->UnitLength_in_cgs;
-  units_time = us->UnitTime_in_cgs;
+  // Set parameter values for chemistry.
+  grackle_data.use_grackle = 1;
+  grackle_data.with_radiative_cooling = 1;
+  /* molecular network with H, He, D
+   From Cloudy table */
+  grackle_data.primordial_chemistry = 0;
+  grackle_data.metal_cooling = 1;  // metal cooling on
+  grackle_data.UVbackground = cooling->uv_background;
+  grackle_data.grackle_data_file = cooling->cloudy_table;
+
+  /* Initialize the chemistry object.
+     a_value is not the true initial a
+     This should get set before any computation */
+  double a_value = 1.;
+
+  if (initialize_chemistry_data(&cooling->units, a_value) == 0) {
+    error("Error in initialize_chemistry_data.");
+  }
 
 #ifdef SWIFT_DEBUG_CHECKS
+  if (GRACKLE_NPART != 1)
+    error("Grackle with multiple particles not implemented");
   float threshold = cooling->GrackleHSShieldingDensityThreshold;
 
   threshold /= phys_const->const_proton_mass;
@@ -200,21 +241,11 @@ static INLINE void cooling_init_backend(
   message("***************************************");
   message("initializing grackle cooling function");
   message("");
-  message("CloudyTable                        = %s",
-          cooling->GrackleCloudyTable);
-  message("UVbackground                       = %d", UVbackground);
-  message("GrackleRedshift                    = %g", cooling->GrackleRedshift);
-  message("GrackleHSShieldingDensityThreshold = %g atom/cm3", threshold);
-#endif
+  cooling_print_backend(cooling);
+  message("Density Self Shielding = %g atom/cm3", threshold);
 
-  if (wrap_init_cooling(cooling->GrackleCloudyTable, UVbackground,
-                        units_density, units_length, units_time,
-                        grackle_chemistry) != 1) {
-    error("Error in initialize_chemistry_data.");
-  }
 
-#ifdef SWIFT_DEBUG_CHECKS
-  grackle_print_data();
+  //grackle_print_data();
   message("");
   message("***************************************");
 #endif
@@ -229,12 +260,90 @@ static INLINE void cooling_print_backend(
     const struct cooling_function_data* cooling) {
 
   message("Cooling function is 'Grackle'.");
-  message("CloudyTable                        = %s",
-          cooling->GrackleCloudyTable);
-  message("UVbackground                       = %d", cooling->UVbackground);
-  message("GrackleRedshift                    = %g", cooling->GrackleRedshift);
-  message("GrackleHSShieldingDensityThreshold = %g atom/cm3",
-          cooling->GrackleHSShieldingDensityThreshold);
+  message("CloudyTable             = %s",
+          cooling->cloudy_table);
+  message("UVbackground            = %d", cooling->uv_background);
+  message("Redshift                = %g", cooling->redshift);
+  message("Density Self Shielding  = %g",
+          cooling->density_self_shielding);
+  message("Units:");
+  message("\tComoving     = %g", cooling->units.comoving_coordinates)
+  message("\tLength  = %g", cooling->units.length_units);
+  message("\tDensity = %g", cooling->units.density_units);
+  message("\tTime         = %g", cooling->units.time_units);
+  message("\tScale Factor = %g", cooling->units.a_units);  
+}
+
+
+/**
+ * @brief print data in grackle struct
+ *
+ * Should only be used for debugging
+ */
+void grackle_print_data() {
+  message("Grackle Data:");
+  message("\t Data file: %s", grackle_data.grackle_data_file);
+  message("\t With grackle: %i", grackle_data.use_grackle);
+  message("\t With radiative cooling: %i", grackle_data.with_radiative_cooling);
+  message("\t With UV background: %i", grackle_data.UVbackground);
+  message("\t With primordial chemistry: %i",
+          grackle_data.primordial_chemistry);
+  message("\t Number temperature bins: %i",
+          grackle_data.NumberOfTemperatureBins);
+  message("\t T = (%g, ..., %g)", grackle_data.TemperatureStart,
+          grackle_data.TemperatureEnd);
+
+  message("Primordial Cloudy");
+  cloudy_print_data(grackle_data.cloudy_primordial, 1);
+  if (grackle_data.metal_cooling) {
+    message("Metal Cooling");
+    cloudy_print_data(grackle_data.cloudy_metal, 0);
+  }
+
+  message("\t Gamma: %g", grackle_data.Gamma);
+
+  /* UVB */
+  if (grackle_data.UVbackground && grackle_data.primordial_chemistry != 0) {
+    struct UVBtable uvb = grackle_data.UVbackground_table;
+    long long N = uvb.Nz;
+    message("\t UV Background");
+    message("\t\t Redshift from %g to %g with %lli steps", uvb.zmin, uvb.zmax,
+            N);
+    message("\t\t z = (%g, ..., %g)", uvb.z[0], uvb.z[N - 1]);
+  }
+}
+
+/**
+ * @brief print data in cloudy struct
+ *
+ * Should only be used for debugging
+ */
+void cloudy_print_data(const cloudy_data c, const int print_mmw) {
+  long long N = c.data_size;
+  message("\t Data size: %lli", N);
+  message("\t Grid rank: %lli", c.grid_rank);
+
+  char msg[200] = "\t Dimension: (";
+  for (long long i = 0; i < c.grid_rank; i++) {
+    char tmp[200] = "%lli%s";
+    if (i == c.grid_rank - 1)
+      sprintf(tmp, tmp, c.grid_dimension[i], ")");
+    else
+      sprintf(tmp, tmp, c.grid_dimension[i], ", ");
+
+    strcat(msg, tmp);
+  }
+  message("%s", msg);
+
+  if (c.heating_data)
+    message("\t Heating: (%g, ..., %g)", c.heating_data[0],
+            c.heating_data[N - 1]);
+  if (c.cooling_data)
+    message("\t Cooling: (%g, ..., %g)", c.cooling_data[0],
+            c.cooling_data[N - 1]);
+  if (c.mmw_data && print_mmw)
+    message("\t Mean molecular weigth: (%g, ..., %g)", c.mmw_data[0],
+            c.mmw_data[N - 1]);
 }
 
 #endif /* SWIFT_COOLING_GRACKLE_H */

src/cooling/grackle/cooling_struct.h

@@ -19,6 +19,9 @@
 #ifndef SWIFT_COOLING_STRUCT_NONE_H
 #define SWIFT_COOLING_STRUCT_NONE_H
 
+/* include grackle */
+#include <grackle.h>
+
 /**
  * @file src/cooling/none/cooling_struct.h
  * @brief Empty infrastructure for the cases without cooling function
@@ -30,16 +33,20 @@
 struct cooling_function_data {
 
   /* Filename of the Cloudy Table */
-  char GrackleCloudyTable[200];
+  char cloudy_table[200];
 
   /* Enable/Disable UV backgroud */
-  int UVbackground;
+  int uv_background;
 
   /* Redshift to use for the UV backgroud (-1 to use cosmological one) */
-  double GrackleRedshift;
+  double redshift;
 
   /* Density Threshold for the shielding */
-  double GrackleHSShieldingDensityThreshold;
+  double density_self_shielding;
+
+  /* unit system */
+  code_units units;
+
 };
 
 /**

src/cooling/grackle/grackle_wrapper.c

@@ -82,18 +82,6 @@ int wrap_init_cooling(char *CloudyTable, int UVbackground, double udensity,
   return 1;
 }
 
-int wrap_set_UVbackground_on() {
-  // The UV background rates is enabled
-  grackle_data.UVbackground = 1;
-  return 1;
-}
-
-int wrap_set_UVbackground_off() {
-  // The UV background rates is disabled
-  grackle_data.UVbackground = 0;
-  return 1;
-}
-
 int wrap_get_cooling_time(double rho, double u, double Z, double a_now,
                           double *coolingtime) {
   gr_float den_factor = 1.0;
@@ -165,63 +153,3 @@ int wrap_do_cooling(double rho, double *u, double dt, double Z, double a_now) {
   return 1;
 }
 
-void grackle_print_data() {
-  message("Grackle Data:");
-  message("\t Data file: %s", grackle_data.grackle_data_file);
-  message("\t With grackle: %i", grackle_data.use_grackle);
-  message("\t With radiative cooling: %i", grackle_data.with_radiative_cooling);
-  message("\t With UV background: %i", grackle_data.UVbackground);
-  message("\t With primordial chemistry: %i",
-          grackle_data.primordial_chemistry);
-  message("\t Number temperature bins: %i",
-          grackle_data.NumberOfTemperatureBins);
-  message("\t T = (%g, ..., %g)", grackle_data.TemperatureStart,
-          grackle_data.TemperatureEnd);
-
-  message("Primordial Cloudy");
-  cloudy_print_data(grackle_data.cloudy_primordial, 1);
-  if (grackle_data.metal_cooling) {
-    message("Metal Cooling");
-    cloudy_print_data(grackle_data.cloudy_metal, 0);
-  }
-
-  message("\t Gamma: %g", grackle_data.Gamma);
-
-  /* UVB */
-  if (grackle_data.UVbackground && grackle_data.primordial_chemistry != 0) {
-    struct UVBtable uvb = grackle_data.UVbackground_table;
-    long long N = uvb.Nz;
-    message("\t UV Background");
-    message("\t\t Redshift from %g to %g with %lli steps", uvb.zmin, uvb.zmax,
-            N);
-    message("\t\t z = (%g, ..., %g)", uvb.z[0], uvb.z[N - 1]);
-  }
-}
-
-void cloudy_print_data(const cloudy_data c, const int print_mmw) {
-  long long N = c.data_size;
-  message("\t Data size: %lli", N);
-  message("\t Grid rank: %lli", c.grid_rank);
-
-  char msg[200] = "\t Dimension: (";
-  for (long long i = 0; i < c.grid_rank; i++) {
-    char tmp[200] = "%lli%s";
-    if (i == c.grid_rank - 1)
-      sprintf(tmp, tmp, c.grid_dimension[i], ")");
-    else
-      sprintf(tmp, tmp, c.grid_dimension[i], ", ");
-
-    strcat(msg, tmp);
-  }
-  message("%s", msg);
-
-  if (c.heating_data)
-    message("\t Heating: (%g, ..., %g)", c.heating_data[0],
-            c.heating_data[N - 1]);
-  if (c.cooling_data)
-    message("\t Cooling: (%g, ..., %g)", c.cooling_data[0],
-            c.cooling_data[N - 1]);
-  if (c.mmw_data && print_mmw)
-    message("\t Mean molecular weigth: (%g, ..., %g)", c.mmw_data[0],
-            c.mmw_data[N - 1]);
-}