[newchem-cpp] lookup_cool_rate1d cleanup part 2

src/clib/initialize_chemistry_data.cpp

@@ -49,18 +49,21 @@ static void show_version(FILE *fp)
   fprintf (fp, "\n");
 }
 
-/**
- * Initialize an empty #gr_interp_grid
- */
-static void initialize_empty_interp_grid_(gr_interp_grid* grid)
-{
-  grid->props.rank = 0;
+/// initialize an empty #gr_interp_grid_props
+static void init_empty_interp_grid_props_(gr_interp_grid_props* props) {
+  props->rank = 0;
   for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++){
-    grid->props.dimension[i] = 0;
-    grid->props.parameters[i] = NULL;
-    grid->props.parameter_spacing[i] = 0.0;
+    props->dimension[i] = 0;
+    props->parameters[i] = nullptr;
+    props->parameter_spacing[i] = 0.0;
   }
-  grid->props.data_size = 0;
+  props->data_size = 0;
+}
+
+/// Initialize an empty #gr_interp_grid
+static void initialize_empty_interp_grid_(gr_interp_grid* grid)
+{
+  init_empty_interp_grid_props_(&(grid->props));
   grid->data=NULL;
 }
 
@@ -397,6 +400,8 @@ extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry,
   my_rates->opaque_storage->kcol_rate_tables = nullptr;
   my_rates->opaque_storage->used_kcol_rate_indices = nullptr;
   my_rates->opaque_storage->n_kcol_rate_indices = 0;
+  init_empty_interp_grid_props_(
+    &my_rates->opaque_storage->h2dust_grain_interp_props);
 
   double co_length_units, co_density_units;
   if (my_units->comoving_coordinates == TRUE) {
@@ -629,6 +634,7 @@ extern "C" int local_free_chemistry_data(chemistry_data *my_chemistry,
     delete my_rates->opaque_storage->kcol_rate_tables;
   }
   delete[] my_rates->opaque_storage->used_kcol_rate_indices;
+  free_interp_grid_props_(&my_rates->opaque_storage->h2dust_grain_interp_props);
   delete my_rates->opaque_storage;
   my_rates->opaque_storage = nullptr;
 

src/clib/initialize_rates.cpp

@@ -250,6 +250,58 @@ int add_h2dust_S_reaction_rate(double **rate_ptr, double units, chemistry_data *
     return GR_SUCCESS;
 }
 
+/// sets up the species-specific h2dust grain coefficient grids
+int setup_h2dust_grain_rates(chemistry_data* my_chemistry,
+                             chemistry_data_storage *my_rates,
+                             double kUnit) {
+
+  //H2 formation on dust grains with C and S compositions
+  if (
+    (add_h2dust_C_reaction_rate(&my_rates->h2dustC, kUnit, my_chemistry)
+      != GR_SUCCESS) ||
+    (add_h2dust_S_reaction_rate(&my_rates->h2dustS, kUnit, my_chemistry)
+      != GR_SUCCESS)
+  ) {
+    return GR_FAIL;
+  }
+
+  // initialize my_rates->opaque_storage->h2dust_grain_interp_props
+  long long n_Tdust = (long long)(my_chemistry->NumberOfDustTemperatureBins);
+  long long n_Tgas = (long long)(my_chemistry->NumberOfTemperatureBins);
+  double* d_Td = (double*)malloc(n_Tdust * sizeof(double));
+  double* d_Tg = (double*)malloc(n_Tgas * sizeof(double));
+
+  const double logtem_start = std::log(my_chemistry->TemperatureStart);
+  const double dlogtem = (std::log(my_chemistry->TemperatureEnd) -
+                          std::log(my_chemistry->TemperatureStart)) /
+                         (double)(my_chemistry->NumberOfTemperatureBins - 1);
+
+  const double logTdust_start = std::log(my_chemistry->DustTemperatureStart);
+  const double dlogTdust = (std::log(my_chemistry->DustTemperatureEnd) -
+                            std::log(my_chemistry->DustTemperatureStart)) /
+                           (double)(my_chemistry->NumberOfDustTemperatureBins - 1);
+
+  for (long long idx = 0; idx < n_Tdust; idx++) {
+    d_Td[idx] = logTdust_start + (double)idx * dlogTdust;
+  }
+  for (long long idx = 0; idx < n_Tgas; idx++) {
+    d_Tg[idx] = logtem_start + (double)idx * dlogtem;
+  }
+
+  gr_interp_grid_props* interp_props
+    = &(my_rates->opaque_storage->h2dust_grain_interp_props);
+  interp_props->rank = 2ll;
+  interp_props->dimension[0] = n_Tdust;
+  interp_props->dimension[1] = n_Tgas;
+  interp_props->parameters[0] = d_Td;
+  interp_props->parameters[1] = d_Tg;
+  interp_props->parameter_spacing[0] = dlogTdust;
+  interp_props->parameter_spacing[0] = dlogtem;
+  interp_props->data_size = n_Tdust*n_Tgas;
+
+  return GR_SUCCESS;
+}
+
 // Down below we define functionality to initialize the table of ordinary
 // collisional rates. If we more fully embraced C++ (and used templates rather
 // than C-style function pointers), this could all be a lot more concise.
@@ -562,9 +614,8 @@ int grackle::impl::initialize_rates(
         //(Equation 9, Wolfire et al., 1995)
         add_reaction_rate(&my_rates->regr, regr_rate, coolingUnits, my_chemistry);
 
-        //H2 formation on dust grains with C and S compositions
-        add_h2dust_C_reaction_rate(&my_rates->h2dustC, kUnit, my_chemistry);
-        add_h2dust_S_reaction_rate(&my_rates->h2dustS, kUnit, my_chemistry);
+        // set up the species-specific h2dust grain coefficient grids
+        setup_h2dust_grain_rates(my_chemistry, my_rates, kUnit);
 
         //Heating of dust by interstellar radiation field, with an arbitrary grain size distribution
         add_scalar_reaction_rate(&my_rates->gamma_isrf2, gamma_isrf2_rate, coolingUnits, my_chemistry);

src/clib/interp_table_utils.h

@@ -21,13 +21,18 @@
 extern "C" {
 #endif
 
+/// Free memory associated with a #gr_interp_grid_props instance
+static inline void free_interp_grid_props_(gr_interp_grid_props* props)
+{
+  for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) {
+    GRACKLE_FREE(props->parameters[i]);
+  }
+}
 
 /// Free memory associated with a #gr_interp_grid
 static inline void free_interp_grid_(gr_interp_grid* grid)
 {
-  for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) {
-    GRACKLE_FREE(grid->props.parameters[i]);
-  }
+  free_interp_grid_props_(&(grid->props));
   GRACKLE_FREE(grid->data);
 }