Moved do_cooling to cooling.h

src/cooling/grackle/cooling.h

@@ -37,7 +37,7 @@
 #include "physical_constants.h"
 #include "units.h"
 
-/* need to rework code if changed */
+/* need to rework (and check) code if changed */
 #define GRACKLE_NPART 1
 
 
@@ -83,26 +83,47 @@ __attribute__((always_inline)) INLINE static double cooling_rate(
     const struct cooling_function_data* restrict cooling,
     struct part* restrict p, float dt) {
 
-  if (cooling->GrackleRedshift == -1) error("TODO time dependant redshift");
+  /* set current time */
+  float scale_factor;
+  if (cooling->redshift == -1)
+    error("TODO time dependant redshift");
+  else
+    scale_factor = 1. / (1. + cooling->redshift);
 
   /* Get current internal energy (dt=0) */
-  double u_old = hydro_get_internal_energy(p);
-  double u_new = u_old;
+  const double energy_before = hydro_get_internal_energy(p);
+  
   /* Get current density */
   const float rho = hydro_get_density(p);
-  /* Actual scaling fractor */
-  const float a_now = 1. / (1. + cooling->GrackleRedshift);
-
+  
   /* 0.02041 (= 1 Zsun in Grackle v2.0, but = 1.5761 Zsun in
      Grackle v2.1) */
-  double Z = 0.02041;
-
-  if (wrap_do_cooling(rho, &u_new, dt, Z, a_now) == 0) {
-    error("Error in do_cooling.\n");
-    return 0;
+  const double Z = 0.02041;
+
+  /* create grackle struct */
+  /* velocities */
+  gr_float x_velocity[GRACKLE_NPART] = {0.0};
+  gr_float y_velocity[GRACKLE_NPART] = {0.0};
+  gr_float z_velocity[GRACKLE_NPART] = {0.0};
+
+  /* particle data */
+  gr_float density[GRACKLE_NPART] = {rho};
+  gr_float metal_density[GRACKLE_NPART] = {Z * density[0]};
+  gr_float energy[GRACKLE_NPART] = {energy_before};
+
+  /* dimensions */
+  int grid_dimension[3] = {GRACKLE_NPART, 0, 0};
+  int grid_start[3] = {0, 0, 0};
+  int grid_end[3] = {0, 0, 0};
+
+  /* solve chemistry with table */
+  if (solve_chemistry_table(&cooling->units, scale_factor, dt, grid_rank, grid_dimension,
+                            grid_start, grid_end, density, energy, x_velocity,
+                            y_velocity, z_velocity, metal_density) == 0) {
+    error("Error in solve_chemistry.");
   }
 
-  return (u_new - u_old) / dt;
+  return (energy[0] - energy_before) / dt;
 }
 
 /**

src/cooling/grackle/grackle_wrapper.c

@@ -31,57 +31,6 @@ gr_float HDI_density[FIELD_SIZE];
 // grid_start and grid_end are used to ignore ghost zones.
 const int grid_rank = 1;
 
-int wrap_init_cooling(char *CloudyTable, int UVbackground, double udensity,
-                      double ulength, double utime, int grackle_chemistry) {
-
-#ifdef SWIFT_DEBUG_CHECKS
-  grackle_verbose = 1;
-#endif
-  double velocity_units;
-
-  // First, set up the units system.
-  // These are conversions from code units to cgs.
-
-  my_units.a_units = 1.0;  // units for the expansion factor (1/1+zi)
-
-  my_units.comoving_coordinates =
-      0; /* so, according to the doc, we assume here all physical quantities to
-            be in proper coordinate (not comobile)  */
-  my_units.density_units = udensity;
-  my_units.length_units = ulength;
-  my_units.time_units = utime;
-  velocity_units =
-      my_units.a_units * my_units.length_units / my_units.time_units;
-  my_units.velocity_units = velocity_units;
-
-  // Second, create a chemistry object for parameters and rate data.
-  if (set_default_chemistry_parameters() == 0) {
-    error("Error in set_default_chemistry_parameters.");
-  }
-
-  // Set parameter values for chemistry.
-  grackle_data.use_grackle = 1;
-  grackle_data.with_radiative_cooling = 1;
-  grackle_data.primordial_chemistry =
-      grackle_chemistry;           // molecular network with H, He, D
-  grackle_data.metal_cooling = 1;  // metal cooling on
-  grackle_data.UVbackground = UVbackground;
-  grackle_data.grackle_data_file = CloudyTable;
-
-  // Finally, initialize the chemistry object.
-  // snl: a_value is not the true initial ae!! This should get set during
-  // update_UVbackground
-  double initial_redshift = 0.;
-  double a_value = 1. / (1. + initial_redshift);
-
-  // Finally, initialize the chemistry object.
-  if (initialize_chemistry_data(&my_units, a_value) == 0) {
-    error("Error in initialize_chemistry_data.");
-  }
-
-  return 1;
-}
-
 int wrap_get_cooling_time(double rho, double u, double Z, double a_now,
                           double *coolingtime) {
   gr_float den_factor = 1.0;