WIP: Refactor Output_float32

src/global/global.cpp

@@ -141,20 +141,6 @@ bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameter
   /* Copy into correct entry in parameters struct */
   if (strcmp(name, "nfile") == 0) {
     parms->nfile = atoi(value);
-  } else if (strcmp(name, "out_float32_density") == 0) {
-    parms->out_float32_density = atoi(value);
-  } else if (strcmp(name, "out_float32_momentum_x") == 0) {
-    parms->out_float32_momentum_x = atoi(value);
-  } else if (strcmp(name, "out_float32_momentum_y") == 0) {
-    parms->out_float32_momentum_y = atoi(value);
-  } else if (strcmp(name, "out_float32_momentum_z") == 0) {
-    parms->out_float32_momentum_z = atoi(value);
-  } else if (strcmp(name, "out_float32_Energy") == 0) {
-    parms->out_float32_Energy = atoi(value);
-#ifdef DE
-  } else if (strcmp(name, "out_float32_GasEnergy") == 0) {
-    parms->out_float32_GasEnergy = atoi(value);
-#endif  // DE
   } else if (strcmp(name, "output_always") == 0) {
     parms->output_always = atoi(value);
 #ifdef MHD

src/global/global.h

@@ -191,23 +191,15 @@ struct Parameters {
   Real gamma;
   char init[MAXLEN];
   int nfile;
-  int n_hydro                = 1;
-  int n_particle             = 1;
-  int n_projection           = 1;
-  int n_rotated_projection   = 1;
-  int n_slice                = 1;
-  int n_out_float32          = 0;
-  int out_float32_density    = 0;
-  int out_float32_momentum_x = 0;
-  int out_float32_momentum_y = 0;
-  int out_float32_momentum_z = 0;
-  int out_float32_Energy     = 0;
-#ifdef DE
-  int out_float32_GasEnergy = 0;
-#endif
-  bool output_always      = false;
-  bool legacy_flat_outdir = false;
-  int n_steps_limit       = -1;  // Note that negative values indicate that there is no limit
+  int n_hydro              = 1;
+  int n_particle           = 1;
+  int n_projection         = 1;
+  int n_rotated_projection = 1;
+  int n_slice              = 1;
+  int n_out_float32        = 0;
+  bool output_always       = false;
+  bool legacy_flat_outdir  = false;
+  int n_steps_limit        = -1;  // Note that negative values indicate that there is no limit
 #ifdef STATIC_GRAV
   int custom_grav = 0;  // flag to set specific static gravity field
 #endif

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;