Introduce FieldInfo

src/grid/field_info.cpp

@@ -0,0 +1,113 @@
+/*! \file
+ *  Define machinery for accessing field information
+ */
+
+#include "field_info.h"
+
+#include <string>
+#include <vector>
+
+#include "../utils/FrozenKeyIdxBiMap.h"
+#include "grid_enum.h"
+
+namespace
+{  // stuff in an anonymous namespace is local to this file
+
+struct PropPack {
+  const char* name;
+  field::Kind kind;
+  field::IOBuf io_buf;
+};
+
+}  // anonymous namespace
+
+/*! list of all field names
+ *
+ *  This must remain synchronized with grid_enum.h
+ */
+static constexpr PropPack pack_arr_[] = {
+    {"density", field::Kind::HYDRO, field::IOBuf::DEVICE},
+    {"momentum_x", field::Kind::HYDRO, field::IOBuf::DEVICE},
+    {"momentum_y", field::Kind::HYDRO, field::IOBuf::DEVICE},
+    {"momentum_z", field::Kind::HYDRO, field::IOBuf::DEVICE},
+    {"Energy", field::Kind::HYDRO, field::IOBuf::DEVICE},
+
+#ifdef SCALAR
+  #ifdef BASIC_SCALAR
+    // we use the name "scalar0" for better consistency with the name recorded during IO
+    {"scalar0", field::Kind::PASSIVE_SCALAR, field::IOBuf::DEVICE},
+  #endif
+
+  #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU)
+    {"HI_density", field::Kind::PASSIVE_SCALAR, field::IOBuf::HOST},
+    {"HII_density", field::Kind::PASSIVE_SCALAR, field::IOBuf::HOST},
+    {"HeI_density", field::Kind::PASSIVE_SCALAR, field::IOBuf::HOST},
+    {"HeII_density", field::Kind::PASSIVE_SCALAR, field::IOBuf::HOST},
+    {"HeIII_density", field::Kind::PASSIVE_SCALAR, field::IOBuf::HOST},
+    {"e_density", field::Kind::PASSIVE_SCALAR, field::IOBuf::HOST},
+    #ifdef GRACKLE_METALS
+    {"metal_density", field::Kind::PASSIVE_SCALAR, field::IOBuf::HOST},
+    #endif
+  #endif
+
+  #ifdef DUST
+    {"dust_density", field::Kind::PASSIVE_SCALAR, field::IOBuf::DEVICE},
+  #endif  // DUST
+
+#endif  // SCALAR
+
+#ifdef MHD
+    {"magnetic_x", field::Kind::MAGNETIC, field::IOBuf::DEVICE},
+    {"magnetic_y", field::Kind::MAGNETIC, field::IOBuf::DEVICE},
+    {"magnetic_z", field::Kind::MAGNETIC, field::IOBuf::DEVICE},
+#endif
+#ifdef DE
+    {"GasEnergy", field::Kind::HYDRO, field::IOBuf::DEVICE}
+#endif
+};
+
+static constexpr int n_fields_ = static_cast<int>(sizeof(pack_arr_) / sizeof(PropPack));
+
+static_assert(n_fields_ == grid_enum::num_fields, "pack_arr_ and grid_enum::num_fields are no longer synchronized");
+
+FieldInfo FieldInfo::create()
+{
+  FieldInfo out;
+
+  // convert n_fields to a vector of std::string
+  std::vector<std::string> v;
+  v.reserve(n_fields_);
+  for (std::size_t i = 0; i < n_fields_; i++) {
+    v.push_back(std::string(pack_arr_[i].name));
+    out.io_buf_.push_back(pack_arr_[i].io_buf);
+    switch (pack_arr_[i].kind) {
+      case field::Kind::HYDRO:
+        out.hydro_field_ids_.push_back(i);
+        break;
+      case field::Kind::PASSIVE_SCALAR:
+        out.scalar_field_ids_.push_back(i);
+        break;
+      case field::Kind::MAGNETIC:
+        out.magnetic_field_ids_.push_back(i);
+        break;
+      default:
+        CHOLLA_ERROR("This branch should be unreachable");
+    }
+  }
+  out.name_id_bimap_ = utils::FrozenKeyIdxBiMap(v);
+  return out;
+}
+
+const std::vector<int>& FieldInfo::get_kind_ids_(field::Kind kind) const
+{
+  switch (kind) {
+    case field::Kind::HYDRO:
+      return hydro_field_ids_;
+    case field::Kind::PASSIVE_SCALAR:
+      return scalar_field_ids_;
+    case field::Kind::MAGNETIC:
+      return magnetic_field_ids_;
+    default:
+      CHOLLA_ERROR("This branch should be unreachable");
+  }
+}

src/grid/field_info.h

@@ -0,0 +1,150 @@
+/*! \file
+ *  Define machinery for accessing field information
+ */
+
+#pragma once
+
+#include <optional>
+#include <string>
+#include <string_view>
+#include <vector>
+
+#include "../utils/FrozenKeyIdxBiMap.h"
+
+namespace field
+{
+
+// note: HYDRO includes GasEnergy (if present)
+enum class Kind { HYDRO, PASSIVE_SCALAR, MAGNETIC };
+
+/*! Specifies which buffer to use for IO */
+enum class IOBuf { HOST, DEVICE };
+
+/*! Specifies centering */
+
+/*! This is a "range" in the C++ 20 sense
+ *
+ *  See \ref FieldInfo::get_id_range for an example
+ */
+class IdRange
+{
+  const std::vector<int>& id_vec_;
+
+ public:
+  explicit IdRange(const std::vector<int>& id_vec) : id_vec_(id_vec) {}
+
+  // the fact that the iterator aliases a const iterator of a std::vector is an
+  // implementation detail
+  using iterator = std::vector<int>::const_iterator;
+
+  iterator begin() const { return id_vec_.begin(); }
+  iterator end() const { return id_vec_.end(); }
+};
+
+}  // namespace field
+
+/*! Dynamically describes the available fields and associated properties
+ */
+class FieldInfo
+{
+  utils::FrozenKeyIdxBiMap name_id_bimap_;
+  std::vector<int> hydro_field_ids_;
+  std::vector<int> scalar_field_ids_;
+  std::vector<int> magnetic_field_ids_;
+  std::vector<field::IOBuf> io_buf_;
+
+  // We make the default-constructor private to force the use of the factory method
+  FieldInfo() = default;
+
+  /*! return a reference to the internal vector of field ids corresponding to
+   *  @ref field::Kind
+   */
+  const std::vector<int>& get_kind_ids_(field::Kind kind) const;
+
+ public:
+  /*! Factory method
+   *
+   *  Ideally, we would make it possible to customize the active scalars, but that's a
+   *  topic for the future.
+   */
+  static FieldInfo create();
+
+  FieldInfo(FieldInfo&&)            = default;
+  FieldInfo& operator=(FieldInfo&&) = default;
+
+  // we delete copy constructor and copy-assignment to prevent accidental copies
+  // (of course move constructors/move assignment remain possible)
+  // In the unlikely event we decide to support copies, this can always change later...
+  FieldInfo(const FieldInfo&)            = delete;
+  FieldInfo& operator=(const FieldInfo&) = delete;
+
+  /*! Get the underlying mapping object between field names and ids */
+  const utils::FrozenKeyIdxBiMap& get_field_id_map() const { return name_id_bimap_; }
+
+  /*! try to lookup the field_id associated with the field_name */
+  std::optional<int> field_id(const char* field_name) const { return name_id_bimap_.find(field_name); }
+  std::optional<int> field_id(std::string_view field_name) const { return name_id_bimap_.find(field_name); }
+
+  /*! try to look up the field name from the field id */
+  std::optional<std::string> field_name(int field_id) const
+  {
+    bool bad_id = (field_id < 0 || field_id >= n_fields());
+    return bad_id ? std::nullopt : std::optional<std::string>{name_id_bimap_.inverse_find(field_id)};
+  }
+
+  /*! try to look up whether the field id refers to a cell-centered field
+   *
+   *  \note
+   *  It may be more useful to return a value that directly specifies whether a field is
+   *  cell-center, x-face-centered, y-face-centered, z-face-centered. It might be
+   *  convenient to specify this with a @ref hydro_utilities::VectorXYZ<int>. For
+   *  example, `{0,0,0}` could represent a cell-centered value and `{1,0,0}`, or maybe
+   *  `{-1, 0, 0}` (we need to think about conventions), could denote a field centered
+   *  on x-faces.
+   */
+  std::optional<bool> is_cell_centered(int field_id)
+  {
+    if (field_id < 0 || field_id >= n_fields()) {
+      return std::nullopt;
+    }
+    for (int id : magnetic_field_ids_) {
+      if (field_id == id) {
+        return std::optional<bool>{false};
+      }
+    }
+    return std::optional<bool>{true};
+  }
+
+  /*! try to look up the IOBuf value associated with a field
+   *
+   *  \note We may want to revisit whether this actually should be tracked by FieldInfo in the future.
+   */
+  std::optional<field::IOBuf> io_buf(int field_id) const
+  {
+    bool bad_id = (field_id < 0 || field_id >= n_fields());
+    return bad_id ? std::nullopt : std::optional<field::IOBuf>{io_buf_[field_id]};
+  }
+
+  /*! Returns the number of fields */
+  int n_fields() const { return static_cast<int>(name_id_bimap_.size()); }
+
+  /*! Returns the number of fields of a given category */
+  int n_fields(field::Kind kind) const { return static_cast<int>(get_kind_ids_(kind).size()); }
+
+  /*! Returns the first field_id corresponding to a passive scalar (if there are any) */
+  std::optional<int> scalar_start() const
+  {
+    return scalar_field_ids_.empty() ? std::nullopt : std::optional<int>{scalar_field_ids_[0]};
+  }
+
+  /*! This returns a "range" over all ids
+   *
+   *  This might be used in a case like the following:
+   *  \code{c++}
+   *  for (int field_id: field_info.get_id_range(field::Kind::HYDRO)) {
+   *    // ...
+   *  }
+   *  \endcode
+   */
+  field::IdRange get_id_range(field::Kind kind) const { return field::IdRange(get_kind_ids_(kind)); }
+};

src/grid/grid3D.cpp

@@ -7,6 +7,7 @@
   #include <hdf5.h>
 #endif
 #include "../global/global.h"
+#include "../grid/field_info.h"
 #include "../grid/grid3D.h"
 #include "../grid/grid_enum.h"       // provides grid_enum
 #include "../hydro/average_cells.h"  // provides Average_Slow_Cells and SlowCellConditionChecker
@@ -39,7 +40,7 @@
 
 /*! \fn Grid3D(void)
  *  \brief Constructor for the Grid. */
-Grid3D::Grid3D(void)
+Grid3D::Grid3D(void) : field_info(FieldInfo::create())
 {
   // set initialization flag to 0
   flag_init = 0;
@@ -112,24 +113,7 @@ Real Grid3D::Calc_Inverse_Timestep()
  *  \brief Initialize the grid. */
 void Grid3D::Initialize(struct Parameters *P)
 {
-  // number of fields to track (default 5 is # of conserved variables)
-  H.n_fields = 5;
-
-// if including passive scalars increase the number of fields
-#ifdef SCALAR
-  H.n_fields += NSCALARS;
-#endif
-
-// if including magnetic fields increase the number of fields
-#ifdef MHD
-  H.n_fields += 3;
-#endif  // MHD
-
-// if using dual energy formalism must track internal energy - always the last
-// field!
-#ifdef DE
-  H.n_fields++;
-#endif
+  H.n_fields = field_info.n_fields();
 
   int nx_in = P->nx;
   int ny_in = P->ny;

src/grid/grid3D.h

@@ -14,6 +14,7 @@
 
 #include "../global/global.h"
 #include "../global/global_cuda.h"
+#include "../grid/field_info.h"
 #include "../io/FnameTemplate.h"
 
 #ifdef HDF5
@@ -50,6 +51,12 @@
   #include "../analysis/analysis.h"
 #endif
 
+// forward declare the DatasetSpec struct
+namespace io
+{
+struct DatasetSpec;
+}  // namespace io
+
 struct Rotation {
   /*! \var nx
    *   \brief Number of pixels in x-dir of rotated, projected image*/
@@ -297,6 +304,9 @@ class Grid3D
    *  \brief Rotation struct for data projections */
   struct Rotation R;
 
+  /*! Describes the mapping between field names and field indices */
+  FieldInfo field_info;
+
 #ifdef GRAVITY
   // Object that contains data for gravity
   Grav3D Grav;
@@ -499,7 +509,7 @@ class Grid3D
 
   /*! \fn void Write_Grid_HDF5(hid_t file_id)
    *  \brief Write the grid to a file, at the current simulation time. */
-  void Write_Grid_HDF5(hid_t file_id);
+  void Write_Grid_HDF5(hid_t file_id, const io::DatasetSpec &h5_dataset_spec);
 
   /*! \fn void Write_Projection_HDF5(hid_t file_id)
    *  \brief Write projected density and temperature data to a file. */

src/grid/grid_enum.h

@@ -22,6 +22,9 @@
 
 // The ": int" forces underlying type to be int.
 
+// if you update this enum, make sure that you also update the pack_arr_ variable from
+// field_info.cpp
+
 namespace grid_enum
 {
 enum : int {