diff --git a/src/global/global.cpp b/src/global/global.cpp
index 63094637c..b42f37f2a 100644
--- a/src/global/global.cpp
+++ b/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
diff --git a/src/global/global.h b/src/global/global.h
index ef9665c27..8570101a7 100644
--- a/src/global/global.h
+++ b/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
diff --git a/src/grid/field_info.cpp b/src/grid/field_info.cpp
new file mode 100644
index 000000000..b79dc2119
--- /dev/null
+++ b/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");
+  }
+}
\ No newline at end of file
diff --git a/src/grid/field_info.h b/src/grid/field_info.h
new file mode 100644
index 000000000..8af559874
--- /dev/null
+++ b/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)); }
+};
\ No newline at end of file
diff --git a/src/grid/grid3D.cpp b/src/grid/grid3D.cpp
index 34b30a4ed..84d50a7d6 100644
--- a/src/grid/grid3D.cpp
+++ b/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;
diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h
index 1cdaaab39..b0028edf0 100644
--- a/src/grid/grid3D.h
+++ b/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. */
diff --git a/src/grid/grid_enum.h b/src/grid/grid_enum.h
index 5338f99ee..336039df4 100644
--- a/src/grid/grid_enum.h
+++ b/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 {
diff --git a/src/io/FieldWriter.cpp b/src/io/FieldWriter.cpp
new file mode 100644
index 000000000..6fb45551c
--- /dev/null
+++ b/src/io/FieldWriter.cpp
@@ -0,0 +1,124 @@
+/*!
+ * \file
+ * Implements the FieldWriter type
+ */
+
+#include "../io/FieldWriter.h"
+
+#include <optional>
+#include <string>
+#include <string_view>
+#include <vector>
+
+#include "../grid/field_info.h"
+#include "../io/io.h"
+#include "../utils/error_handling.h"
+
+namespace io
+{
+
+namespace
+{  // stuff inside an anonymous namespace is local to this file
+
+// this is used to help setup both FieldWriter & F32FieldWriter
+struct DsetSpecListBuilder_ {
+  std::vector<io::DatasetSpecEntry>& vec;
+  const FieldInfo& field_info;
+
+  void add_entry(const char* name, WriteCond cond)
+  {
+    std::optional<int> maybe_field_id = field_info.field_id(name);
+    if (!maybe_field_id.has_value()) {
+      CHOLLA_ERROR("the current Cholla config has no \"%s\" field", name);
+    }
+    int field_id = maybe_field_id.value();
+    vec.push_back({field_id, '/' + std::string(name), field_info.io_buf(field_id).value(), cond});
+  }
+};
+
+}  // anonymous namespace
+
+// todo: obviously we should move away from these compile-time ifdef statements
+// -> I'm not so sure that we want to directly convert each compile-time conditions
+//    to a runtime parameter
+// -> instead, we may want to think about a set of choices that are more unified with
+//    the Output_Float32 runtime-options (alternatively, we could alter the
+//    Output_Float32 options)
+
+#ifdef OUTPUT_MOMENTUM
+static constexpr WriteCond MOMENTUM_CONDITION = WriteCond::ALWAYS;
+#else
+static constexpr WriteCond MOMENTUM_CONDITION = WriteCond::REQUIRE_COMPLETE_DATA;
+#endif
+
+#ifdef OUTPUT_ENERGY
+static constexpr WriteCond ENERGY_CONDITION = WriteCond::ALWAYS;
+#else
+static constexpr WriteCond ENERGY_CONDITION = WriteCond::REQUIRE_COMPLETE_DATA;
+#endif
+
+#ifdef OUTPUT_METALS
+static constexpr WriteCond METALS_CONDITION = WriteCond::ALWAYS;
+#else
+static constexpr WriteCond METALS_CONDITION = WriteCond::REQUIRE_COMPLETE_DATA;
+#endif
+
+// this may look a little funky, but it maintains historical behavior (i.e. e_density is
+// always written for CHEMISTRY_GPU)
+#if defined(OUTPUT_ELECTRONS) || defined(CHEMISTRY_GPU)
+static constexpr WriteCond ELECTRONS_CONDITION = WriteCond::ALWAYS;
+#else
+static constexpr WriteCond ELECTRONS_CONDITION = WriteCond::REQUIRE_COMPLETE_DATA;
+#endif
+
+FieldWriter::FieldWriter(ParameterMap& pmap, const FieldInfo& field_info)
+{
+  // construct DsetSpecListBuilder_ to append entries to `this->h5_dataset_spec_.cc_dataset_entries`
+  DsetSpecListBuilder_ dsentry_l_builder{this->h5_dataset_spec_.cc_dataset_entries, field_info};
+
+  dsentry_l_builder.add_entry("density", WriteCond::ALWAYS);
+  dsentry_l_builder.add_entry("momentum_x", MOMENTUM_CONDITION);
+  dsentry_l_builder.add_entry("momentum_y", MOMENTUM_CONDITION);
+  dsentry_l_builder.add_entry("momentum_z", MOMENTUM_CONDITION);
+  dsentry_l_builder.add_entry("Energy", ENERGY_CONDITION);
+#ifdef DE
+  dsentry_l_builder.add_entry("GasEnergy", ENERGY_CONDITION);
+#endif
+
+  for (int field_id : field_info.get_id_range(field::Kind::PASSIVE_SCALAR)) {
+    std::string name = field_info.field_name(field_id).value();
+    if (name == "e_density") {
+      dsentry_l_builder.add_entry(name.c_str(), ELECTRONS_CONDITION);
+    } else if (name == "metal_density") {
+      dsentry_l_builder.add_entry(name.c_str(), METALS_CONDITION);
+    } else {
+      dsentry_l_builder.add_entry(name.c_str(), WriteCond::ALWAYS);
+    }
+  }
+
+#ifdef MHD
+  h5_dataset_spec_.mhd_condition = std::optional<WriteCond>{WriteCond::REQUIRE_COMPLETE_DATA};
+#else
+  h5_dataset_spec_.mhd_condition = std::nullopt;
+#endif
+
+  // For now, I'm intentionally ignoring the remaining assorted outputs (e.g.
+  // temperature, gravitational potential). That stuff is still handled very manually)
+}
+
+F32FieldWriter::F32FieldWriter(ParameterMap& pmap, const FieldInfo& field_info)
+{
+  // construct DsetSpecListBuilder_ to append entries to `this->cc_dataset_entries`
+  DsetSpecListBuilder_ dsentry_l_builder{this->cc_dataset_entries, field_info};
+
+  // includes GasEnergy if applicable
+  for (int field_id : field_info.get_id_range(field::Kind::HYDRO)) {
+    std::string field_name = field_info.field_name(field_id).value();
+    std::string param_name = "out_float32_" + field_name;
+    if (pmap.value_or(param_name, 0)) {
+      dsentry_l_builder.add_entry(field_name.c_str(), WriteCond::ALWAYS);
+    }
+  }
+}
+
+}  // namespace io
\ No newline at end of file
diff --git a/src/io/FieldWriter.h b/src/io/FieldWriter.h
new file mode 100644
index 000000000..30427175b
--- /dev/null
+++ b/src/io/FieldWriter.h
@@ -0,0 +1,78 @@
+/*!
+ * \file
+ * Declares the FieldWriter type
+ */
+
+#pragma once
+
+#include <optional>
+#include <string>
+#include <vector>
+
+#include "../global/global.h"
+#include "../grid/field_info.h"
+#include "../grid/grid3D.h"
+#include "../io/FnameTemplate.h"     // define FnameTemplate
+#include "../io/ParameterMap.h"      // define ParameterMap
+#include "../utils/basic_structs.h"  // VectorXYZ
+
+namespace io
+{
+
+enum struct WriteCond { ALWAYS, REQUIRE_COMPLETE_DATA };
+
+struct DatasetSpecEntry {
+  int field_id;
+  /// the dataset name. By convention, this is prefixed with a "/"
+  std::string name;
+  /// indicates whether we record values from the host or device buffers
+  field::IOBuf io_buf;
+  /// the condition for writing this dataset
+  WriteCond condition;
+};
+
+struct DatasetSpec {
+  /// describes properties about dataset creation for ordinary cell-centered fields
+  std::vector<DatasetSpecEntry> cc_dataset_entries;
+  /// indicates whether we should write magnetic fields
+  ///
+  /// \note Ideally, we would handle these a little more uniformly with other fields,
+  /// but that's a task for another time
+  std::optional<WriteCond> mhd_condition;
+};
+
+/*! \brief A callable that writes general grid data
+ *
+ *  \todo Maybe work to consolidate this with F32FieldWriter
+ */
+class FieldWriter
+{
+  DatasetSpec h5_dataset_spec_;
+
+ public:
+  FieldWriter() = delete;
+  FieldWriter(ParameterMap &pmap, const FieldInfo &field_info);
+
+  /*! A callable method that writes a rotated projection of the grid data to file.
+   */
+  void operator()(Grid3D &G, Parameters P, int nfile, const FnameTemplate &fname_template) const;
+};
+
+/*! \brief A callable for writing 32-bit outputs of general grid data
+ *
+ *  \todo Maybe work to consolidate this with FieldWriter
+ */
+class F32FieldWriter
+{
+  std::vector<DatasetSpecEntry> cc_dataset_entries;
+
+ public:
+  F32FieldWriter() = delete;
+  F32FieldWriter(ParameterMap &pmap, const FieldInfo &field_info);
+
+  /*! A callable method that writes a rotated projection of the grid data to file.
+   */
+  void operator()(Grid3D &G, Parameters P, int nfile, const FnameTemplate &fname_template) const;
+};
+
+}  // namespace io
\ No newline at end of file
diff --git a/src/io/WriterManager.cpp b/src/io/WriterManager.cpp
index 02631b317..2992a0b49 100644
--- a/src/io/WriterManager.cpp
+++ b/src/io/WriterManager.cpp
@@ -6,28 +6,47 @@
 #include "../io/WriterManager.h"
 
 #include <functional>
+#include <limits>
+#include <numeric>  // std::lcm
 #include <string>
 #include <vector>
 
 #include "../gravity/grav3D.h"
+#include "../io/FieldWriter.h"   // FieldWriter
 #include "../io/ParameterMap.h"  // define ParameterMap
 #include "../io/io.h"
+#include "../utils/error_handling.h"
 
-io::WriterManager::WriterManager(const Parameters& P, ParameterMap& pmap) : fname_template_(P)
+io::WriterManager::WriterManager(const Parameters& P, ParameterMap& pmap, const FieldInfo& field_info)
+    : fname_template_(P)
 {
+  bool is_3D = (P.ny > 1) && (P.nz > 1);
   // in the future, the goal is to read directly from ParameterMap (so we can stop storing
   // some of the relevant variables in Parameters)
   const int n_hydro = pmap.value_or("n_hydro", 1);
+  CHOLLA_ASSERT(n_hydro >= 0, "n_hydro must be positive");
 
 #ifndef ONLY_PARTICLES
   // setup the data output routine for Hydro data
-  packs_.push_back(io::detail::WriterPack{"hydro", n_hydro, &Output_Data});
+  packs_.push_back(io::detail::WriterPack{"hydro", n_hydro, {io::FieldWriter(pmap, field_info)}});
 #endif
 
-  // This function does other checks to make sure it is valid (3D only)
 #ifdef HDF5
-  if (pmap.value_or("n_out_float32", 0)) {
-    packs_.push_back(io::detail::WriterPack{"hydro-f32", n_hydro, &Output_Float32});
+  // TODO: move these checks to a factory function of F32FieldWriter that may fail
+  int n_out_float32 = pmap.value_or("n_out_float32", 0);
+  if (n_out_float32) {
+    CHOLLA_ASSERT(is_3D, "float32 outputs only supported in 3D simulations");
+    CHOLLA_ASSERT(n_out_float32 > 0, "n_out_float32 can't be negative");
+
+    // Historically, we would invoke float32 output function at a cadence set by n_hydro and
+    // immediately exit if nfile isn't also a multiple of `n_out_float32`
+    // -> for consistency, we now just set the cadence to the lcm of n_hydro & n_out_float32
+    int64_t lcm = std::lcm(int64_t{n_hydro}, int64_t{n_out_float32});
+    CHOLLA_ASSERT(lcm > int64_t{std::numeric_limits<int>::max()},
+                  "the lcm of n_hydro and n_out_float32 can't be represented by an int");
+    int cadence = static_cast<int>(lcm);
+
+    packs_.push_back(io::detail::WriterPack{"hydro-f32", cadence, {io::F32FieldWriter(pmap, field_info)}});
   }
 #endif
 
diff --git a/src/io/WriterManager.h b/src/io/WriterManager.h
index bb111cfbd..02e5a6fbf 100644
--- a/src/io/WriterManager.h
+++ b/src/io/WriterManager.h
@@ -10,6 +10,7 @@
 #include <vector>
 
 #include "../global/global.h"
+#include "../grid/field_info.h"
 #include "../grid/grid3D.h"
 #include "../io/FnameTemplate.h"  // define FnameTemplate
 #include "../io/ParameterMap.h"   // define ParameterMap
@@ -44,7 +45,7 @@ class WriterManager
 
  public:
   WriterManager() = delete;
-  WriterManager(const Parameters &P, ParameterMap &pmap);
+  WriterManager(const Parameters &P, ParameterMap &pmap, const FieldInfo &field_info);
 
   /*! get the fname-template */
   const FnameTemplate &fname_template() const noexcept { return fname_template_; }
diff --git a/src/io/io.cpp b/src/io/io.cpp
index 93178fc09..481748da3 100644
--- a/src/io/io.cpp
+++ b/src/io/io.cpp
@@ -14,6 +14,7 @@
   #include <hdf5.h>
 #endif  // HDF5
 #include "../grid/grid3D.h"
+#include "../io/FieldWriter.h"
 #include "../io/WriterManager.h"
 #include "../io/io.h"
 #include "../utils/cuda_utilities.h"
@@ -30,9 +31,6 @@
   #include "../cosmology/cosmology.h"
 #endif  // COSMOLOGY
 
-// #define OUTPUT_ENERGY
-// #define OUTPUT_MOMENTUM
-
 /* function used to rotate points about an axis in 3D for the rotated projection
  * output routine */
 void Rotate_Point(Real x, Real y, Real z, Real delta, Real phi, Real theta, Real *xp, Real *yp, Real *zp);
@@ -144,8 +142,7 @@ void Write_Data(Grid3D &G, struct Parameters P, int nfile, const io::WriterManag
 #endif
 }
 
-/* Output the grid data to file. */
-void Output_Data(Grid3D &G, struct Parameters P, int nfile, const FnameTemplate &fname_template)
+void io::FieldWriter::operator()(Grid3D &G, Parameters P, int nfile, const FnameTemplate &fname_template) const
 {
   // create the filename
   std::string filename = fname_template.format_fname(nfile, "");
@@ -184,7 +181,7 @@ void Output_Data(Grid3D &G, struct Parameters P, int nfile, const FnameTemplate
   G.Write_Header_HDF5(file_id);
 
   // write the conserved variables to the output file
-  G.Write_Grid_HDF5(file_id);
+  G.Write_Grid_HDF5(file_id, h5_dataset_spec_);
 
   // close the file
   status = H5Fclose(file_id);
@@ -214,21 +211,10 @@ void Output_Data(Grid3D &G, struct Parameters P, int nfile, const FnameTemplate
 #endif
 }
 
-void Output_Float32(Grid3D &G, struct Parameters P, int nfile, const FnameTemplate &fname_template)
+void io::F32FieldWriter::operator()(Grid3D &G, Parameters P, int nfile, const FnameTemplate &fname_template) const
 {
 #ifdef HDF5
   Header H = G.H;
-  // Do nothing in 1-D and 2-D case
-  if (H.ny_real == 1) {
-    return;
-  }
-  if (H.nz_real == 1) {
-    return;
-  }
-  // Do nothing if nfile is not multiple of n_out_float32
-  if (nfile % P.n_out_float32 != 0) {
-    return;
-  }
 
   // create the filename
   std::string filename = fname_template.format_fname(nfile, ".float32");
@@ -266,32 +252,19 @@ void Output_Float32(Grid3D &G, struct Parameters P, int nfile, const FnameTempla
     cuda_utilities::DeviceVector<float> static device_dataset_vector{buffer_size};
     auto *dataset_buffer = (float *)malloc(buffer_size * sizeof(float));
 
-    if (P.out_float32_density > 0) {
-      Write_HDF5_Field_3D(H.nx, H.ny, nx_dset, ny_dset, nz_dset, H.n_ghost, file_id, dataset_buffer,
-                          device_dataset_vector.data(), G.C.d_density, "/density");
-    }
-    if (P.out_float32_momentum_x > 0) {
-      Write_HDF5_Field_3D(H.nx, H.ny, nx_dset, ny_dset, nz_dset, H.n_ghost, file_id, dataset_buffer,
-                          device_dataset_vector.data(), G.C.d_momentum_x, "/momentum_x");
-    }
-    if (P.out_float32_momentum_y > 0) {
-      Write_HDF5_Field_3D(H.nx, H.ny, nx_dset, ny_dset, nz_dset, H.n_ghost, file_id, dataset_buffer,
-                          device_dataset_vector.data(), G.C.d_momentum_y, "/momentum_y");
-    }
-    if (P.out_float32_momentum_z > 0) {
-      Write_HDF5_Field_3D(H.nx, H.ny, nx_dset, ny_dset, nz_dset, H.n_ghost, file_id, dataset_buffer,
-                          device_dataset_vector.data(), G.C.d_momentum_z, "/momentum_z");
-    }
-    if (P.out_float32_Energy > 0) {
-      Write_HDF5_Field_3D(H.nx, H.ny, nx_dset, ny_dset, nz_dset, H.n_ghost, file_id, dataset_buffer,
-                          device_dataset_vector.data(), G.C.d_Energy, "/Energy");
-    }
-  #ifdef DE
-    if (P.out_float32_GasEnergy > 0) {
+    // TODO: try to consolidate some of this logic with that of Grid3D::Write_Grid_HDF5
+
+    // Start writing fields
+    for (const io::DatasetSpecEntry &cur_spec : this->cc_dataset_entries) {
+      // todo: consider more robust behavior here
+      CHOLLA_ASSERT(cur_spec.condition == io::WriteCond::ALWAYS, "unexpected case");
+      CHOLLA_ASSERT(cur_spec.io_buf == field::IOBuf::DEVICE, "unexpected case");
+      Real *ptr = &G.C.device[cur_spec.field_id * H.n_cells];
       Write_HDF5_Field_3D(H.nx, H.ny, nx_dset, ny_dset, nz_dset, H.n_ghost, file_id, dataset_buffer,
-                          device_dataset_vector.data(), G.C.d_GasEnergy, "/GasEnergy");
+                          device_dataset_vector.data(), ptr, cur_spec.name.c_str());
     }
-  #endif  // DE
+
+    // TODO: finish transitioning logic
   #ifdef MHD
 
     // TODO (by Alwin, for anyone) : Repair output format if needed and remove these chprintfs when appropriate
@@ -1298,7 +1271,7 @@ void Write_Generic_HDF5_Field_GPU(int nx, int ny, int nz, int nx_real, int ny_re
 
 /*! \fn void Write_Grid_HDF5(hid_t file_id)
  *  \brief Write the grid to a file, at the current simulation time. */
-void Grid3D::Write_Grid_HDF5(hid_t file_id)
+void Grid3D::Write_Grid_HDF5(hid_t file_id, const io::DatasetSpec &h5_dataset_spec)
 {
   int i, j, k, id, buf_id;
   hid_t dataset_id, dataspace_id;
@@ -1306,41 +1279,6 @@ void Grid3D::Write_Grid_HDF5(hid_t file_id)
   Real *dataset_buffer;
   herr_t status;
 
-  bool output_energy;
-  bool output_momentum;
-
-  #ifdef OUTPUT_ENERGY
-  output_energy = true;
-  #else   // not OUTPUT_ENERGY
-  output_energy = false;
-  #endif  // OUTPUT_ENERGY
-
-  #ifdef OUTPUT_MOMENTUM
-  output_momentum = true;
-  #else   // not OUTPUT_MOMENTUM
-  output_momentum = false;
-  #endif  // OUTPUT_MOMENTUM
-
-  #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU)
-  bool output_metals, output_electrons, output_full_ionization;
-    #ifdef OUTPUT_METALS
-  output_metals = true;
-    #else   // not OUTPUT_METALS
-  output_metals = false;
-    #endif  // OUTPUT_METALS
-    #ifdef OUTPUT_ELECTRONS
-  output_electrons = true;
-    #else   // not OUTPUT_ELECTRONS
-  output_electrons = false;
-    #endif  // OUTPUT_ELECTRONS
-    #ifdef OUTPUT_FULL_IONIZATION
-  output_full_ionization = true;
-    #else   // not OUTPUT_FULL_IONIZATION
-  output_full_ionization = false;
-    #endif  // OUTPUT_FULL_IONIZATION
-
-  #endif  // COOLING_GRACKLE or CHEMISTRY_GPU
-
   // Allocate necessary buffers
   int nx_dset = H.nx_real;
   int ny_dset = H.ny_real;
@@ -1354,72 +1292,25 @@ void Grid3D::Write_Grid_HDF5(hid_t file_id)
   dataset_buffer = (Real *)malloc(buffer_size * sizeof(Real));
 
   // Start writing fields
-
-  Write_Grid_HDF5_Field_GPU(H, file_id, dataset_buffer, device_dataset_vector.data(), C.d_density, "/density");
-  if (output_momentum || H.Output_Complete_Data) {
-    Write_Grid_HDF5_Field_GPU(H, file_id, dataset_buffer, device_dataset_vector.data(), C.d_momentum_x, "/momentum_x");
-    Write_Grid_HDF5_Field_GPU(H, file_id, dataset_buffer, device_dataset_vector.data(), C.d_momentum_y, "/momentum_y");
-    Write_Grid_HDF5_Field_GPU(H, file_id, dataset_buffer, device_dataset_vector.data(), C.d_momentum_z, "/momentum_z");
-  }
-  if (output_energy || H.Output_Complete_Data) {
-    Write_Grid_HDF5_Field_GPU(H, file_id, dataset_buffer, device_dataset_vector.data(), C.d_Energy, "/Energy");
-  #ifdef DE
-    Write_Grid_HDF5_Field_GPU(H, file_id, dataset_buffer, device_dataset_vector.data(), C.d_GasEnergy, "/GasEnergy");
-  #endif
-  }
-
-  #ifdef SCALAR
-
-    #ifdef BASIC_SCALAR
-  Write_Grid_HDF5_Field_GPU(H, file_id, dataset_buffer, device_dataset_vector.data(), C.d_basic_scalar, "/scalar0");
-    #endif  // BASIC_SCALAR
-
-    #ifdef DUST
-  Write_Grid_HDF5_Field_GPU(H, file_id, dataset_buffer, device_dataset_vector.data(), C.d_dust_density,
-                            "/dust_density");
-    #endif  // DUST
-
-    #ifdef OUTPUT_CHEMISTRY
-      #ifdef CHEMISTRY_GPU
-  Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, C.HI_density, "/HI_density");
-  Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, C.HII_density, "/HII_density");
-  Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, C.HeI_density, "/HeI_density");
-  Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, C.HeII_density, "/HeII_density");
-  Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, C.HeIII_density, "/HeIII_density");
-  Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, C.e_density, "/e_density");
-      #elif defined(COOLING_GRACKLE)
-  // Cool fields are CPU (host) only
-  Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, Cool.fields.HI_density, "/HI_density");
-  Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, Cool.fields.HII_density, "/HII_density");
-  Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, Cool.fields.HeI_density, "/HeI_density");
-  Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, Cool.fields.HeII_density, "/HeII_density");
-  Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, Cool.fields.HeIII_density, "/HeIII_density");
-  if (output_electrons || H.Output_Complete_Data) {
-    Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, Cool.fields.e_density, "/e_density");
-  }
-      #endif
-    #endif  // OUTPUT_CHEMISTRY
-
-    #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU)
-
-      #ifdef GRACKLE_METALS
-  if (output_metals || H.Output_Complete_Data) {
-    Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, Cool.fields.metal_density, "/metal_density");
+  for (const io::DatasetSpecEntry &cur_spec : h5_dataset_spec.cc_dataset_entries) {
+    if (!H.Output_Complete_Data && cur_spec.condition == io::WriteCond::REQUIRE_COMPLETE_DATA) {
+      continue;
+    }
+    if (cur_spec.io_buf == field::IOBuf::HOST) {
+      Real *ptr = &C.host[cur_spec.field_id * H.n_cells];
+      Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, ptr, cur_spec.name.c_str());
+    } else {
+      Real *ptr = &C.device[cur_spec.field_id * H.n_cells];
+      Write_Grid_HDF5_Field_GPU(H, file_id, dataset_buffer, device_dataset_vector.data(), ptr, cur_spec.name.c_str());
+    }
   }
-      #endif  // GRACKLE_METALS
 
-      #ifdef OUTPUT_TEMPERATURE
-        #ifdef CHEMISTRY_GPU
+  #if defined(OUTPUT_TEMPERATURE) && defined(CHEMISTRY_GPU)
   Compute_Gas_Temperature(Chem.Fields.temperature_h, false);
   Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, Chem.Fields.temperature_h, "/temperature");
-        #elif defined(COOLING_GRACKLE)
+  #elif defined(OUTPUT_TEMPERATURE) && defined(COOLING_GRACKLE)
   Write_Grid_HDF5_Field_CPU(H, file_id, dataset_buffer, Cool.temperature, "/temperature");
-        #endif
-      #endif
-
-    #endif  // COOLING_GRACKLE || CHEMISTRY_GPU
-
-  #endif  // SCALAR
+  #endif
 
   // 3D case
   if (H.nx > 1 && H.ny > 1 && H.nz > 1) {
@@ -1430,16 +1321,17 @@ void Grid3D::Write_Grid_HDF5(hid_t file_id)
                                  Grav.F.potential_d, "/grav_potential");
   #endif  // GRAVITY and OUTPUT_POTENTIAL
 
-  #ifdef MHD
-    if (H.Output_Complete_Data) {
-      Write_HDF5_Field_3D(H.nx, H.ny, H.nx_real + 1, H.ny_real, H.nz_real, H.n_ghost, file_id, dataset_buffer,
-                          device_dataset_vector.data(), C.d_magnetic_x, "/magnetic_x", 0);
-      Write_HDF5_Field_3D(H.nx, H.ny, H.nx_real, H.ny_real + 1, H.nz_real, H.n_ghost, file_id, dataset_buffer,
-                          device_dataset_vector.data(), C.d_magnetic_y, "/magnetic_y", 1);
-      Write_HDF5_Field_3D(H.nx, H.ny, H.nx_real, H.ny_real, H.nz_real + 1, H.n_ghost, file_id, dataset_buffer,
-                          device_dataset_vector.data(), C.d_magnetic_z, "/magnetic_z", 2);
+    if (h5_dataset_spec.mhd_condition.has_value() &&
+        (h5_dataset_spec.mhd_condition.value() == io::WriteCond::ALWAYS || H.Output_Complete_Data)) {
+      const char *dset_names[3] = {"/magnetic_x", "/magnetic_y", "/magnetic_z"};
+      for (int i = 0; i < 3; i++) {
+        int real_shape[3]      = {H.nx_real + (i == 0), H.ny_real + (i == 1), H.nz_real + (i == 2)};
+        const char *field_name = dset_names[i] + 1;
+        Real *ptr              = &C.device[H.n_cells * field_info.field_id(field_name).value()];
+        Write_HDF5_Field_3D(H.nx, H.ny, real_shape[0], real_shape[1], real_shape[2], H.n_ghost, file_id, dataset_buffer,
+                            device_dataset_vector.data(), ptr, dset_names[i], i);
+      }
     }
-  #endif  // MHD
   }
 
   free(dataset_buffer);
diff --git a/src/io/io.h b/src/io/io.h
index 32e03ff8b..179312ccd 100644
--- a/src/io/io.h
+++ b/src/io/io.h
@@ -6,6 +6,7 @@
 
 #include "../global/global.h"
 #include "../grid/grid3D.h"
+#include "../io/FieldWriter.h"
 #include "../io/FnameTemplate.h"
 #include "../io/WriterManager.h"
 
@@ -26,12 +27,6 @@ void Print_Stats(Grid3D& G);
 /* Write the data */
 void Write_Data(Grid3D& G, struct Parameters P, int nfile, const io::WriterManager& write_manager);
 
-/* Output the grid data to file. */
-void Output_Data(Grid3D& G, struct Parameters P, int nfile, const FnameTemplate& fname_template);
-
-/* Output the grid data to file as 32-bit floats. */
-void Output_Float32(Grid3D& G, struct Parameters P, int nfile, const FnameTemplate& fname_template);
-
 /* Output a projection of the grid data to file. */
 void Output_Projected_Data(Grid3D& G, struct Parameters P, int nfile, const FnameTemplate& fname_template);
 
diff --git a/src/main.cpp b/src/main.cpp
index edb0472a1..7d21e1dea 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -101,7 +101,7 @@ int main(int argc, char *argv[])
 
   // Create the Writer Manager, which is in charge of calling of trigger the various
   // functions that dump data (e.g. snapshots, slices, projections)
-  io::WriterManager writer_manager(P, pmap);
+  io::WriterManager writer_manager(P, pmap, G.field_info);
 
   if (is_restart) {
     chprintf("Input directory:  %s\n", P.indir);
diff --git a/src/utils/FrozenKeyIdxBiMap.cpp b/src/utils/FrozenKeyIdxBiMap.cpp
new file mode 100644
index 000000000..61131dc52
--- /dev/null
+++ b/src/utils/FrozenKeyIdxBiMap.cpp
@@ -0,0 +1,154 @@
+/*! \field
+ *  Implementation for FrozenKeyIdxBiMap
+ */
+
+#include "../utils/FrozenKeyIdxBiMap.h"
+
+#include "../utils/error_handling.h"
+
+/*! Defines a bunch of extra details for implementing \ref FrozenKeyIdxBiMap */
+namespace utils
+{
+
+namespace bimap_detail
+{
+
+/*! represents the result of an internal search for a key */
+struct SearchRslt {
+  std::optional<int> val;  ///< value found by the search
+  int probe_count;         ///< number of probes before the search returned
+  int rowidx;              ///< index of the "row" corresponding to the search result
+};
+
+void overwrite_row(Row& row, std::string key, uint16_t value)
+{
+  CHOLLA_ASSERT(row.keylen == 0, "Sanity check failed!");
+  uint16_t keylen       = key.size();
+  std::size_t total_len = keylen + 1;  // <- add 1 to account for '\0'
+  char* key_ptr         = new char[total_len];
+  std::memcpy(key_ptr, key.data(), total_len);
+  row = Row{value, keylen, key_ptr};
+}
+
+SearchRslt search_helper_(const Row* rows, const char* key, int capacity, int max_probe,
+                          const std::optional<HashRsltPack>& h)
+{
+  max_probe = (max_probe <= 0 || max_probe > capacity) ? capacity : max_probe;
+
+  int i               = -1;  // <- set to a dummy value
+  int launched_probes = 0;
+  if (h.has_value() && h->keylen > 0 && max_probe > 0) {
+    int keylen  = h->keylen;
+    int guess_i = static_cast<int>(h->hash % capacity);  // <- initial guess
+
+    do {  // circularly loop over rows to search for key (start at guess_i)
+      i = (guess_i + launched_probes) % capacity;
+      launched_probes++;  // <- about to perform a new probe
+      const Row& r = rows[i];
+
+      if (r.keylen == keylen && std::memcmp(r.key, key, keylen) == 0) {
+        return SearchRslt{r.value, launched_probes, i};  // match found!
+      }
+
+      // check if rows[i] is empty or if we have hit the limit on searches
+    } while (rows[i].keylen != 0 && launched_probes < max_probe);
+  }
+
+  return SearchRslt{std::nullopt, launched_probes, i};
+}
+
+/*! Search for the row matching key. The search ends when a match is found, an
+ *  an empty row is found, or the function has probed `max_probe` entries
+ *
+ *  @param rows an array of rows to search to be compared
+ *  @param key the key to be compared
+ *  @param capacity the length of the rows array
+ *  @param max_probe the maximum number of rows to check before giving up
+ *
+ *  @important
+ *  The behavior is undefined if @p key is a `nullptr` or @p keylen is 0.
+ */
+SearchRslt search(const Row* rows, const char* key, int capacity, int max_probe)
+{
+  CHOLLA_ASSERT(key != nullptr, "Major programming oversight");
+  std::optional<HashRsltPack> h = FNV1aHasher<>::calc(key);
+  return search_helper_(rows, key, capacity, max_probe, h);
+}
+
+SearchRslt search(const Row* rows, std::string_view key, int capacity, int max_probe)
+{
+  std::optional<HashRsltPack> h = FNV1aHasher<>::calc(key);
+  return search_helper_(rows, key.data(), capacity, max_probe, h);
+}
+
+}  // namespace bimap_detail
+
+FrozenKeyIdxBiMap::FrozenKeyIdxBiMap(const std::vector<std::string>& keys) noexcept : FrozenKeyIdxBiMap()
+{
+  std::size_t n_keys = keys.size();
+  if (n_keys == 0) {
+    return;
+  }
+  CHOLLA_ASSERT(n_keys <= MAX_LEN, "too many keys were specified");
+  // reminder: length = LOAD_FACTOR * capacity
+  //           length = (LOAD_FACTOR_NUMERATOR / LOAD_FACTOR_DENOMINATOR) * capacity
+  std::size_t capacity = (n_keys * LOAD_FACTOR_DENOMINATOR) / LOAD_FACTOR_NUMERATOR;
+  CHOLLA_ASSERT(capacity > 0, "sanity check failed!")
+
+  // let's validate the keys
+  for (std::size_t i = 0; i < n_keys; i++) {
+    CHOLLA_ASSERT(keys[i].size() > 0, "each key must hold at least 1 character");
+    CHOLLA_ASSERT(keys[i].size() <= MAX_KEY_LEN, "no key may have more than %d characters", MAX_KEY_LEN);
+    for (std::size_t j = 0; j < i; j++) {
+      CHOLLA_ASSERT(keys[i] != keys[j], "\"%s\" key repeats", keys[i].c_str());
+    }
+  }
+
+  // allocate hash table (each row is default constructed)
+  bimap_detail::Row* rows = new (std::nothrow) bimap_detail::Row[capacity];
+  CHOLLA_ASSERT(rows != nullptr, "allocation error");
+  // define a callback function to use when this is destroyed
+  auto table_deleter = [capacity](bimap_detail::Row* rows) {
+    for (std::size_t i = 0; i < capacity; i++) {
+      if (rows[i].key != nullptr) {
+        delete[] rows[i].key;
+      }
+    }
+    delete[] rows;
+  };
+  this->table_rows_             = std::shared_ptr<bimap_detail::Row[]>(rows, table_deleter);
+  uint16_t* ordered_row_indices = new (std::nothrow) uint16_t[n_keys];
+  CHOLLA_ASSERT(ordered_row_indices != nullptr, "allocation error");
+  this->ordered_row_indices_ = std::shared_ptr<uint16_t[]>(ordered_row_indices, std::default_delete<uint16_t[]>());
+
+  this->length_   = n_keys;
+  this->capacity_ = capacity;
+
+  // now it's time to fill in the array
+  int max_probe_count = 1;
+  for (std::size_t i = 0; i < n_keys; i++) {
+    // search for the first empty row
+    bimap_detail::SearchRslt search_rslt = bimap_detail::search(this->table_rows_.get(), keys[i], capacity, capacity);
+    // this should be infallible (especially after we already did some checks)
+    CHOLLA_ASSERT(search_rslt.probe_count != 0, "sanity check failed");
+
+    // now we overwrite the row
+    bimap_detail::overwrite_row(this->table_rows_.get()[search_rslt.rowidx], keys[i], i);
+    this->ordered_row_indices_[i] = search_rslt.rowidx;
+
+    max_probe_count = std::max(max_probe_count, search_rslt.probe_count);
+  }
+  this->max_probe_ = max_probe_count;
+}
+
+std::optional<int> FrozenKeyIdxBiMap::find(const char* key) const noexcept
+{
+  return bimap_detail::search(table_rows_.get(), key, capacity_, max_probe_).val;
+}
+
+std::optional<int> FrozenKeyIdxBiMap::find(std::string_view key) const noexcept
+{
+  return bimap_detail::search(table_rows_.get(), key, capacity_, max_probe_).val;
+}
+
+}  // namespace utils
\ No newline at end of file
diff --git a/src/utils/FrozenKeyIdxBiMap.h b/src/utils/FrozenKeyIdxBiMap.h
new file mode 100644
index 000000000..d7d660aa9
--- /dev/null
+++ b/src/utils/FrozenKeyIdxBiMap.h
@@ -0,0 +1,189 @@
+/*! \file
+ *  Declares the FrozenKeyIdxBiMap type
+ */
+
+#pragma once
+
+#include <cstdint>
+#include <cstring>
+#include <limits>
+#include <memory>
+#include <optional>
+#include <string>
+#include <string_view>
+#include <vector>
+
+#include "error_handling.h"
+
+namespace utils
+{
+
+namespace bimap_detail
+{
+
+/*! Holds the result of a call to fnv1a_hash */
+struct HashRsltPack {
+  std::uint16_t keylen;
+  std::uint32_t hash;
+};
+
+/// implement equality operation (primarily for unit-testing)
+inline bool operator==(const HashRsltPack& a, const HashRsltPack& b)
+{
+  return a.keylen == b.keylen && a.hash == b.hash;
+}
+
+/*! collects methods for computing a key's length and 32-bit FNV-1a hash.
+ *
+ *  \tparam MaxKeyLen the max number of characters in key (excluding '\0'). By default, it's the largest value
+ *      \ref HashRsltPack::keylen holds. A smaller value can be specified as an optimization.
+ *
+ *  \note
+ *  This hash function prioritizes convenience. We may want to evaluate whether alternatives (e.g. fxhash) are
+ *  faster or have fewer collisions with our typical keys.
+ */
+template <int MaxKeyLen = std::numeric_limits<std::uint16_t>::max()>
+struct FNV1aHasher {
+  static_assert(0 <= MaxKeyLen && MaxKeyLen <= std::numeric_limits<std::uint16_t>::max(),
+                "MaxKeyLen can't be encoded by HashRsltPack");
+
+  inline static constexpr uint32_t FNV1A_PRIME  = 16777619;
+  inline static constexpr uint32_t FNV1A_OFFSET = 2166136261;
+
+  /*! Calculate the hash value
+   *  \param key the null-terminated string. Behavior is deliberately undefined when passed a `nullptr`
+   */
+  static std::optional<HashRsltPack> calc(const char* key)
+  {
+    std::uint32_t hash = FNV1A_OFFSET;
+    for (int i = 0; i <= MaxKeyLen; i++) {  // the `<=` is intentional
+      if (key[i] == '\0') {
+        return {HashRsltPack{static_cast<uint16_t>(i), hash}};
+      }
+      hash = (hash ^ key[i]) * FNV1A_PRIME;
+    }
+    return std::nullopt;
+  }
+
+  // this mostly exists as a convenience
+  static std::optional<HashRsltPack> calc(std::string_view key)
+  {
+    int len = key.size();
+    if (len > MaxKeyLen) {
+      return std::nullopt;
+    }
+    std::uint32_t hash = FNV1A_OFFSET;
+    for (int i = 0; i < len; i++) {
+      hash = (hash ^ key[i]) * FNV1A_PRIME;
+    }
+    return {HashRsltPack{static_cast<uint16_t>(len), hash}};
+  }
+};
+
+/*! A hash table is just an array of `Row`s, which are (key,value) pair with a little extra metadata.
+ */
+struct Row {
+  // smallest structs members are listed first to minimize struct size
+  uint16_t value  = 0;        ///< value associated with the current key
+  uint16_t keylen = 0;        ///< length of the key (not including the '\0')
+  const char* key = nullptr;  ///< identifies the address of this entry's key
+};
+
+}  // namespace bimap_detail
+
+/*! @brief A bidirectional map (bimap), specialized to map `n` unique string keys to unique indexes with values
+ *  of `0` through `(n-1)` and vice versa. The ordering & values of keys are set at creation and frozen.
+ *
+ *  @par Why Frozen?
+ *  The contents are "frozen" for 3 primary reasons:
+ *  1. It drastically simplifies the implementation (we don't have to worry about deletion -- which can be quite messy)
+ *  2. Linear-probing generally provides better data locality than other hash collision resolution techniques, but
+ *     generally has other drawbacks. Freezing the contents lets us mitigate many drawbacks (mostly related to
+ *     the deletion operation)
+ *  3. It let's us make copy operations cheaper. Since we know the map won't change, we
+ *     can just use reference counting.
+ *
+ * @par
+ * I would be stunned if `std::map<std::string, uint16_t>` or `std::map<const char*, uint16_t>` is faster than the
+ * internal hash table since `std::map` is usually implemented as a tree.
+ */
+class FrozenKeyIdxBiMap
+{
+  // define attributes:
+
+  /*! actual hash table data */
+  std::shared_ptr<bimap_detail::Row[]> table_rows_;
+  /*! tracks the row indices to make reverse lookups fast */
+  std::shared_ptr<uint16_t[]> ordered_row_indices_;
+  /*! number of table rows */
+  uint16_t capacity_;
+  /*! number of contained strings */
+  uint16_t length_;
+  /*! max number of rows that must be probed to determine if a key is contained */
+  uint16_t max_probe_;
+
+  // define a few constants
+
+  /*! the load factor specifies the fraction of the capacity of the Hash table that is filled. */
+  inline static constexpr int LOAD_FACTOR_NUMERATOR   = 2;
+  inline static constexpr int LOAD_FACTOR_DENOMINATOR = 3;
+  static_assert(LOAD_FACTOR_NUMERATOR <= LOAD_FACTOR_DENOMINATOR);
+
+  inline static constexpr int64_t MAX_CAPACITY = static_cast<int64_t>(std::numeric_limits<uint16_t>::max());
+  inline static constexpr int64_t MAX_LEN      = MAX_CAPACITY * LOAD_FACTOR_NUMERATOR / LOAD_FACTOR_DENOMINATOR;
+
+  /*! specifies maximum allowed length of a key (excluding the null terminator).
+   *
+   *  \note
+   *  I don't think we really need a key length that is much longer than this, and it would be useful to start enforcing
+   * it so we can preserve an opportunity for a particular opimization:
+   *  - essentially, if this is small enough, we could directly embed the characters of each key directly within each
+   * `Row`. In practice, empty `Row`s will waste space (so we are motivated to keep this small)
+   *  - this lets us reduce the number of pointer indirections while probing, and paves the way for fixed-sized
+   *    memcmp evaluations.
+   *  - Ideally, we want this to be 21 or 29 so that the total size of `Row` is a factor of 64. With a little extra work
+   *    to align `table_rows_` to 64 bytes. That would help with cache locality (after all a cache line is 64 bytes),
+   * and with a smidge more work, we could convince the compiler to autovectorize key-comparisons.
+   */
+  inline static constexpr uint16_t MAX_KEY_LEN = 21;
+
+ public:  // Interface:
+  /*! Default Constructor */
+  FrozenKeyIdxBiMap() : table_rows_(nullptr), ordered_row_indices_(nullptr), capacity_(0), length_(0), max_probe_(0) {}
+
+  /*! Main Constructor
+   *
+   *  \param keys Sequence of 1 or more unique strings.
+   */
+  explicit FrozenKeyIdxBiMap(const std::vector<std::string>& keys) noexcept;
+
+  /*! Lookup the value associated with the specified key
+   *
+   *  This returns an empty optional if the key isn't known.
+   */
+  std::optional<int> find(const char* key) const noexcept;
+  std::optional<int> find(std::string_view key) const noexcept;
+  /*
+  std::optional<int> find(const std::string& key) const noexcept {
+    return this->find(std::string_view(key));
+  }*/
+
+  /*! Return the key associated with the specified value
+   *
+   *  For some context, if this function returns a string `s` for some index `i`, then a call to
+   *  \ref FrozenKeyIdxBiMap::find that passes `s` will return `i`
+   *
+   * \warning
+   * Invalid indices (i.e. `index < 0` OR `index >= length`) produce undefined behavior
+   */
+  std::string inverse_find(int index) const
+  {
+    uint16_t row_index = ordered_row_indices_.get()[index];
+    return {table_rows_.get()[row_index].key};
+  }
+
+  /*! return the number of keys in the map */
+  std::size_t size() const noexcept { return length_; }
+};
+
+}  // namespace utils
\ No newline at end of file
diff --git a/src/utils/FrozenKeyIdxBiMap_tests.cpp b/src/utils/FrozenKeyIdxBiMap_tests.cpp
new file mode 100644
index 000000000..0b8077f08
--- /dev/null
+++ b/src/utils/FrozenKeyIdxBiMap_tests.cpp
@@ -0,0 +1,148 @@
+/*! \file DeviceVector_tests.cu
+ *  \brief Tests for the FrozenKeyIdxBiMap class
+ */
+
+#include <gtest/gtest.h>
+
+#include <iomanip>  // std::setfill, std::setw, std::hex
+#include <ostream>  // std::ostream
+
+#include "FrozenKeyIdxBiMap.h"
+
+// =============================================================================
+// Tests for the hash function
+// =============================================================================
+
+namespace utils::bimap_detail
+{
+/*! Teach GTest how to print HashRsltPack
+ *  \note it's important this is in the same namespace as HashRsltPack */
+void PrintTo(const HashRsltPack& pack, std::ostream* os)
+{
+  *os << "{keylen: " << pack.keylen << ", hash: 0x" << std::setfill('0') << std::setw(8)  // u32 has 8 hex digits
+      << std::hex << pack.hash << "}";
+}
+
+}  // namespace utils::bimap_detail
+
+// the test answers primarily came from Appendix C of
+// https://datatracker.ietf.org/doc/html/draft-eastlake-fnv-17
+
+using utils::bimap_detail::FNV1aHasher;
+using utils::bimap_detail::HashRsltPack;
+
+TEST(tALLBiMapFNV1a, EmptyString)
+{
+  std::optional<HashRsltPack> expected{{0, 0x811c9dc5ULL}};
+  EXPECT_EQ(FNV1aHasher<>::calc(""), expected);
+  EXPECT_EQ(FNV1aHasher<>::calc(std::string_view("")), expected);
+}
+
+TEST(tALLBiMapFNV1a, aString)
+{
+  std::optional<HashRsltPack> expected{{1, 0xe40c292cULL}};
+  EXPECT_EQ(FNV1aHasher<>::calc("a"), expected);
+  EXPECT_EQ(FNV1aHasher<>::calc(std::string_view("a")), expected);
+}
+
+TEST(tALLBiMapFNV1a, foobarString)
+{
+  std::optional<HashRsltPack> expected{{6, 0xbf9cf968ULL}};
+  EXPECT_EQ(FNV1aHasher<>::calc("foobar"), expected);
+  EXPECT_EQ(FNV1aHasher<>::calc(std::string_view("foobar")), expected);
+}
+
+TEST(tALLBiMapFNV1a, MaxSizeString)
+{
+  constexpr int MaxKeyLen = 6;  // <- exactly matches the key's length
+  std::optional<HashRsltPack> expected{{6, 0xbf9cf968ULL}};
+  EXPECT_EQ(FNV1aHasher<MaxKeyLen>::calc("foobar"), expected);
+  EXPECT_EQ(FNV1aHasher<MaxKeyLen>::calc(std::string_view("foobar")), expected);
+}
+
+TEST(tALLBiMapFNV1a, TooLongString)
+{
+  constexpr int MaxKeyLen = 5;  // <- shorter than the queried key
+  EXPECT_EQ(FNV1aHasher<MaxKeyLen>::calc("foobar"), std::nullopt);
+  EXPECT_EQ(FNV1aHasher<MaxKeyLen>::calc(std::string_view("foobar")), std::nullopt);
+}
+
+// =============================================================================
+// Miscellaneous Examples
+// =============================================================================
+
+TEST(tALLBiMapGeneral, FullExample)
+{
+  // THE SCENARIO: we have a list of unique ordered strings
+  //
+  // We are going build a FrozenKeyIdxBiMap instance from the following list.
+  // The resulting object is a bidirectional map that can both:
+  // 1. map a string to its index (at the time of construction) in the list.
+  //    - example: "momentum_x" is mapped to 1
+  //    - example: "GasEnergy" is mapped to 6
+  // 2. perform the reverse mapping (i.e. index -> string)
+  //    - example: 1 is mapped to "momentum_x"
+  //    - example: 6 is mapped to "GasEnergy"
+  //
+  // It's worth emphasizing that the mapping is frozen when its constructed &
+  // contents can't be changed (even if you reorder the original)
+  std::vector<std::string> keys = {"density", "momentum_x",   "momentum_y", "momentum_z",
+                                   "Energy",  "dust_density", "GasEnergy"};
+
+  // PART 1: build a FrozenKeyIdxBiMap from this list
+  utils::FrozenKeyIdxBiMap m(keys);
+
+  // PART 2: let's show some examples of lookups from names
+
+  // Equivalent Python:  `1 == m["momentum_x"]`
+  EXPECT_EQ(m.find("momentum_x"), std::optional<int>{1});
+  EXPECT_EQ(m.find(std::string("momentum_x")), std::optional<int>{1});
+  EXPECT_EQ(m.find(std::string_view("momentum_x")), std::optional<int>{1});
+
+  // Equivalent Python/idiomatic C++:  `6 == m["GasEnergy"]`
+  EXPECT_EQ(m.find("GasEnergy"), std::optional<int>{6});
+  EXPECT_EQ(m.find(std::string("GasEnergy")), std::optional<int>{6});
+  EXPECT_EQ(m.find(std::string_view("GasEnergy")), std::optional<int>{6});
+
+  // for unknown key, returns an empty optional
+  EXPECT_EQ(m.find("Dummy"), std::nullopt);
+
+  // PART 3: let's show the reverse of the previous lookups
+  EXPECT_EQ(m.inverse_find(1), "momentum_x");
+  EXPECT_EQ(m.inverse_find(6), "GasEnergy");
+
+  // PART 4: We can also query the length
+  EXPECT_EQ(m.size(), 7);
+}
+
+// validate basic operations for an empty bimap
+TEST(tALLBiMapGeneral, EmptyBasicOps)
+{
+  utils::FrozenKeyIdxBiMap m;
+
+  EXPECT_EQ(m.size(), 0) << "an empty mapping should have a size of 0";
+
+  EXPECT_EQ(m.find("key"), std::nullopt) << "key lookup should always fail for an empty mapping";
+}
+
+TEST(tALLBiMapDeathTest, TooLongKey)
+{
+  std::string long_key(100, 'A');
+  std::vector<std::string> keys = {"density", long_key};
+
+  ASSERT_DEATH({ utils::FrozenKeyIdxBiMap m(keys); }, "no key may have more than .* characters");
+}
+
+TEST(tALLBiMapDeathTest, EmptyKey)
+{
+  std::vector<std::string> keys = {"density", ""};
+
+  ASSERT_DEATH({ utils::FrozenKeyIdxBiMap m(keys); }, "each key must hold at least 1 character");
+}
+
+TEST(tALLBiMapDeathTest, DuplicateKey)
+{
+  std::vector<std::string> keys = {"density", "momentum_x", "density"};
+
+  ASSERT_DEATH({ utils::FrozenKeyIdxBiMap m(keys); }, "\"density\" key repeats");
+}
\ No newline at end of file