[newchem-cpp] generalize experimental ratequery API

.clang-format-ignore

@@ -64,17 +64,17 @@ src/include/grackle_float.h.in
 src/include/grackle_misc.h
 src/include/grackle_rate_functions.h
 src/include/grackle_types.h
-tests/unit/grtest_cmd.cpp
-tests/unit/grtest_cmd.hpp
-tests/unit/grtest_os.cpp
-tests/unit/grtest_os.hpp
-tests/unit/grtest_utils.cpp
-tests/unit/grtest_utils.hpp
+tests/grtestutils/cmd.cpp
+tests/grtestutils/cmd.hpp
+tests/grtestutils/os.cpp
+tests/grtestutils/os.hpp
+tests/grtestutils/utils.cpp
+tests/grtestutils/utils.hpp
+tests/grtestutils/googletest/check_allclose.hpp
 tests/unit/test_chemistry_struct_synced.cpp
 tests/unit/test_ghost_zone.cpp
 tests/unit/test_interpolators_comparisons.cpp
 tests/unit/test_linalg.cpp
 tests/unit/test_status_reporting.cpp
 tests/unit/test_unit_interpolators_g.cpp
-tests/unit/utest_helpers.hpp
 

src/clib/CMakeLists.txt

@@ -120,10 +120,11 @@ add_library(Grackle_Grackle
   lookup_cool_rates1d.hpp
   make_consistent.cpp make_consistent.hpp
   opaque_storage.hpp
-  rate_utils.cpp
+  ratequery.cpp ratequery.hpp
   solve_chemistry.cpp
   scale_fields.cpp scale_fields.hpp
   solve_rate_cool_g-cpp.cpp solve_rate_cool_g-cpp.h
+  support/SimpleVec.hpp
   step_rate_gauss_seidel.hpp
   step_rate_newton_raphson.hpp
   time_deriv_0d.hpp

src/clib/Make.config.objects

@@ -49,7 +49,7 @@ OBJS_CONFIG_LIB = \
         calc_tdust_3d.lo \
         calc_grain_size_increment_1d.lo \
         rate_functions.lo \
-        rate_utils.lo \
+        ratequery.lo \
         gaussj_g.lo \
         utils.lo \
         utils-cpp.lo \

src/clib/chemistry_solver_funcs.hpp

@@ -1,6 +1,11 @@
-// See LICENSE file for license and copyright information
-
-/// @file chemistry_solver_funcs.hpp
+//===----------------------------------------------------------------------===//
+//
+// See the LICENSE file for license and copyright information
+// SPDX-License-Identifier: NCSA AND BSD-3-Clause
+//
+//===----------------------------------------------------------------------===//
+///
+/// @file
 /// @brief Defines chemistry reaction related functions invoked by the
 ///     grackle solver in order to integrate the species densities over time.
 ///
@@ -14,6 +19,8 @@
 ///     be decoupled from the derivative calculation for primoridial species
 ///   - it may also make sense to further divide logic by the kinds of species
 ///     that are affected (e.g. primordial vs grains)
+///
+//===----------------------------------------------------------------------===//
 
 #ifndef CHEMISTRY_SOLVER_FUNCS_HPP
 #define CHEMISTRY_SOLVER_FUNCS_HPP
@@ -579,7 +586,9 @@ inline void species_density_updates_gauss_seidel(
 
         if ( (my_chemistry->metal_chemistry == 1)  && 
              (itmask_metal[i] != MASK_FALSE) )  {
-          scoef = scoef;
+          // we comment out the following line that assigns scoef to itself since
+          // it has no practical impact and produces a compiler warning
+          // scoef = scoef;
           acoef = acoef
               + kcol_buf.data[CollisionalRxnLUT::kz44][i] *   CII(i,j,k) / 12.
               + kcol_buf.data[CollisionalRxnLUT::kz45][i] *   OII(i,j,k) / 16.
@@ -1932,7 +1941,9 @@ inline void species_density_derivatives_0d(
 
     if ((my_chemistry->metal_chemistry == 1)  && 
         (itmask_metal[0] != MASK_FALSE))  {
-      scoef = scoef;
+      // we comment out the following line that assigns scoef to itself since
+      // it has no practical impact and produces a compiler warning
+      // scoef = scoef;
       acoef = acoef
           + kcr_buf.data[CollisionalRxnLUT::kz44][0]    *   CII        / 12.
           + kcr_buf.data[CollisionalRxnLUT::kz45][0]    *   OII        / 16.

src/clib/initialize_chemistry_data.cpp

@@ -28,6 +28,8 @@
 #include "internal_types.hpp" // drop_CollisionalRxnRateCollection
 #include "opaque_storage.hpp" // gr_opaque_storage
 #include "phys_constants.h"
+#include "ratequery.hpp"
+#include "status_reporting.h"
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -229,9 +231,15 @@ static void initialize_empty_chemistry_data_storage_struct(chemistry_data_storag
   my_rates->opaque_storage = NULL;
 }
 
-extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry,
-                                               chemistry_data_storage *my_rates,
-                                               code_units *my_units)
+/// core logic of local_initialize_chemistry_data_
+///
+/// @note
+/// This has been separated from local_initialize_chemistry_data to ensure that
+/// any memory allocations tracked by reg_builder can be appropriately freed
+/// (this is somewhat unavoidable in C++ without destructors)
+static int local_initialize_chemistry_data_(
+    chemistry_data *my_chemistry, chemistry_data_storage *my_rates,
+    code_units *my_units, grackle::impl::ratequery::RegBuilder* reg_builder)
 {
 
   /* Better safe than sorry: Initialize everything to NULL/0 */
@@ -403,6 +411,7 @@ extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry,
   init_empty_interp_grid_props_(
     &my_rates->opaque_storage->h2dust_grain_interp_props);
   my_rates->opaque_storage->grain_species_info = nullptr;
+  my_rates->opaque_storage->registry = nullptr;
 
   double co_length_units, co_density_units;
   if (my_units->comoving_coordinates == TRUE) {
@@ -418,7 +427,8 @@ extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry,
 
   // Compute rate tables.
   if (grackle::impl::initialize_rates(my_chemistry, my_rates, my_units,
-                                      co_length_units, co_density_units)
+                                      co_length_units, co_density_units,
+                                      reg_builder)
       != GR_SUCCESS) {
     fprintf(stderr, "Error in initialize_rates.\n");
     return GR_FAIL;
@@ -456,6 +466,16 @@ extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry,
   /* store a copy of the initial units */
   my_rates->initial_units = *my_units;
 
+  // initialize the registry
+  if (grackle::impl::ratequery::RegBuilder_misc_recipies(reg_builder,
+                                                         my_chemistry)
+      != GR_SUCCESS){
+    return GrPrintAndReturnErr("error in RegBuilder_misc_recipies");
+  }
+  my_rates->opaque_storage->registry = new grackle::impl::ratequery::Registry(
+    grackle::impl::ratequery::RegBuilder_consume_and_build(reg_builder)
+  );
+
   if (grackle_verbose) {
     time_t timer;
     char tstr[80];
@@ -503,6 +523,20 @@ extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry,
   return GR_SUCCESS;
 }
 
+
+extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry,
+                                               chemistry_data_storage *my_rates,
+                                               code_units *my_units)
+{
+  namespace rate_q = grackle::impl::ratequery;
+  rate_q::RegBuilder reg_builder = rate_q::new_RegBuilder();
+
+  int out = local_initialize_chemistry_data_(my_chemistry, my_rates, my_units,
+                                             &reg_builder);
+  rate_q::drop_RegBuilder(&reg_builder);
+  return out;
+}
+
 extern "C" int initialize_chemistry_data(code_units *my_units)
 {
   if (local_initialize_chemistry_data(grackle_data, &grackle_rates,
@@ -659,6 +693,15 @@ extern "C" int local_free_chemistry_data(chemistry_data *my_chemistry,
     delete my_rates->opaque_storage->grain_species_info;
   }
 
+  if (my_rates->opaque_storage->registry != nullptr) {
+    // delete contents of registry
+    grackle::impl::ratequery::drop_Registry(
+      my_rates->opaque_storage->registry
+    );
+    // delete registry, itself
+    delete my_rates->opaque_storage->registry;
+  }
+
   delete my_rates->opaque_storage;
   my_rates->opaque_storage = nullptr;
 

src/clib/initialize_rates.cpp

@@ -385,12 +385,79 @@ int init_kcol_rate_tables(
   return GR_SUCCESS;
 }
 
+
+/// a function that encodes the algorithm for looking up the `i`th pointer
+/// collisional reaction rate from an instance of `chemistry_data_storage`.
+///
+/// This is intended to be registered with RegBuilder in order to provide
+/// access to the various collisional reaction rates.
+///
+/// @param my_rates The object from which the rate Entry is loaded
+/// @param i the index of the queried rate
+grackle::impl::ratequery::Entry get_CollisionalRxn_Entry
+    (chemistry_data_storage* my_rates, int i) {
+  // sanity check! (this shouldn't actually happen)
+  if ((my_rates == nullptr) || (my_rates->opaque_storage == nullptr) ||
+      (my_rates->opaque_storage->kcol_rate_tables == nullptr)) {
+    return grackle::impl::ratequery::mk_invalid_Entry();
+  }
+
+  // import new_Entry into the current scope (so we don't need the full name)
+  using ::grackle::impl::ratequery::new_Entry;
+
+  double** data = my_rates->opaque_storage->kcol_rate_tables->data;
+
+  // this implementation leverages the following properties of the
+  // CollisionalRxnLUT enum:
+  // - each entry of collisional_rxn_rate_members.def has a corresponding
+  //   enumeration-constant
+  // - the very first enumeration constant has a value of 0 (since a value
+  //   wasn't explicitly specified)
+  // - the value of each other enumeration constants is 1 larger than the value
+  //   of the previous value (if a value isn't explicitly specified)
+  // - CollisionalRxnLUT::NUM_ENTRIES specifies the number of other enumeration
+  //   constants (excluding CollisionalRxnLUT::NUM_ENTRIES) in the enum
+  switch (i) {
+#define TO_STR(s) #s
+#define ENTRY(NAME)                                                            \
+  case CollisionalRxnLUT::NAME: { return new_Entry(data[i], TO_STR(NAME)); }
+#include "collisional_rxn_rate_members.def"
+
+#undef ENTRY
+#undef TO_STR
+    default: {
+      return grackle::impl::ratequery::mk_invalid_Entry();
+    }
+  }
+}
+
+
+/// a function that encodes the algorithm for looking up the k13dd
+/// collisional reaction rate from an instance of `chemistry_data_storage`.
+///
+/// This is intended to be registered with RegBuilder in order to provide
+/// access to the various collisional reaction rates.
+///
+/// @param my_rates The object from which the rate Entry is loaded
+/// @param i the index of the queried rate
+grackle::impl::ratequery::Entry get_k13dd_Entry(
+    chemistry_data_storage* my_rates, int i) {
+  if (i == 0) {
+    double* ptr = (my_rates == nullptr) ? nullptr : my_rates->k13dd;
+    return grackle::impl::ratequery::new_Entry(ptr, "k13dd");
+  } else {
+    return grackle::impl::ratequery::mk_invalid_Entry();
+  }
+}
+
+
 } // anonymous namespace
 
 //Definition of the initialise_rates function.
 int grackle::impl::initialize_rates(
   chemistry_data *my_chemistry, chemistry_data_storage *my_rates,
-  code_units *my_units, double co_length_unit, double co_density_unit)
+  code_units *my_units, double co_length_unit, double co_density_unit,
+  ratequery::RegBuilder* reg_builder)
 { 
     // TODO: we REALLY need to do an error check that
     //    my_chemistry->NumberOfTemperatureBins >= 2
@@ -493,6 +560,15 @@ int grackle::impl::initialize_rates(
         return GR_FAIL;
       }
 
+      // register the recipe for looking up the "standard" collisional rates
+      if (grackle::impl::ratequery::RegBuilder_recipe_1d(
+              reg_builder, CollisionalRxnLUT::NUM_ENTRIES,
+              &get_CollisionalRxn_Entry, my_chemistry->NumberOfTemperatureBins
+          ) != GR_SUCCESS) {
+        return GrPrintAndReturnErr("error registering standard collisional "
+                                   "reaction rates");
+      }
+
         //--------Calculate coefficients for density-dependent collisional H2 dissociation rate--------
         //
         // idt = 0 calculates coefficients for direct collisional dissociation idt = 1 
@@ -508,6 +584,14 @@ int grackle::impl::initialize_rates(
         //          coeff7(idt=0, TbinFinal), coeff1(idt=1, Tbin1), ..., coeff7(idt=1, TbinFinal)}
         add_k13dd_reaction_rate(&my_rates->k13dd, kUnit, my_chemistry);
 
+        // maybe k13dd should be considered multi-dimensional?
+        if (grackle::impl::ratequery::RegBuilder_recipe_1d(
+                reg_builder, 1, &get_k13dd_Entry,
+                my_chemistry->NumberOfTemperatureBins * 14) != GR_SUCCESS) {
+          return GrPrintAndReturnErr("error registering k13dd rate");
+        }
+
+
         //H2 formation on dust grains requires loop over the dust temperature.
         add_h2dust_reaction_rate(&my_rates->h2dust, kUnit, my_chemistry);
 

src/clib/initialize_rates.hpp

@@ -14,12 +14,14 @@
 #define INITIALIZE_RATES_HPP
 
 #include "grackle.h"
+#include "ratequery.hpp"
 
 namespace grackle::impl {
 
 int initialize_rates(chemistry_data* my_chemistry,
                      chemistry_data_storage* my_rates, code_units* my_units,
-                     double co_length_unit, double co_density_unit);
+                     double co_length_unit, double co_density_unit,
+                     ratequery::RegBuilder* reg_builder);
 
 }  // namespace grackle::impl
 

src/clib/opaque_storage.hpp

@@ -16,6 +16,7 @@
 #include "grackle.h"
 #include "dust/grain_species_info.hpp"
 #include "internal_types.hpp"
+#include "ratequery.hpp"
 
 /// a struct that used to wrap some private storage details
 ///
@@ -95,6 +96,9 @@ struct gr_opaque_storage {
   /// > calculations). An alternative would be to briefly initialize an
   /// > instance during setup and then repack the data.
   grackle::impl::GrainSpeciesInfo* grain_species_info;
+
+  /// used to implement the experimental ratequery machinery
+  grackle::impl::ratequery::Registry* registry;
 };
 
 #endif /* OPAQUE_STORAGE_HPP */