From fd6a52923faab8ca1ab07f5c60c52aa20fd938d4 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 20 Jan 2026 10:38:21 -0500 Subject: [PATCH 01/16] A bunch of the initial changes --- src/utils/FrozenKeyIdxBiMap.cpp | 150 ++++++++++++++++++++ src/utils/FrozenKeyIdxBiMap.h | 189 ++++++++++++++++++++++++++ src/utils/FrozenKeyIdxBiMap_tests.cpp | 148 ++++++++++++++++++++ 3 files changed, 487 insertions(+) create mode 100644 src/utils/FrozenKeyIdxBiMap.cpp create mode 100644 src/utils/FrozenKeyIdxBiMap.h create mode 100644 src/utils/FrozenKeyIdxBiMap_tests.cpp diff --git a/src/utils/FrozenKeyIdxBiMap.cpp b/src/utils/FrozenKeyIdxBiMap.cpp new file mode 100644 index 000000000..ddfc709dc --- /dev/null +++ b/src/utils/FrozenKeyIdxBiMap.cpp @@ -0,0 +1,150 @@ +/*! \field + * Implementation for FrozenKeyIdxBiMap + */ + +#include "FrozenKeyIdxBiMap.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 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& 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(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 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 h = FNV1aHasher<>::calc(key); + return search_helper_(rows, key.data(), capacity, max_probe, h); +} + +} // namespace bimap_detail + +FrozenKeyIdxBiMap::FrozenKeyIdxBiMap(const std::vector& 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 bimap_detail::Row[capacity]; + // 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(rows, table_deleter); + uint16_t* ordered_row_indices = new uint16_t[n_keys]; + this->ordered_row_indices_ = std::shared_ptr(ordered_row_indices, std::default_delete()); + + 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 FrozenKeyIdxBiMap::find(const char* key) const noexcept +{ + return bimap_detail::search(table_rows_.get(), key, capacity_, max_probe_).val; +} + +std::optional 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..ac109d504 --- /dev/null +++ b/src/utils/FrozenKeyIdxBiMap.h @@ -0,0 +1,189 @@ +/*! \file + * Declares the FrozenKeyIdxBiMap type + */ + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#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 ::max()> +struct FNV1aHasher { + static_assert(0 <= MaxKeyLen && MaxKeyLen <= std::numeric_limits::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 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(i), hash}}; + } + hash = (hash ^ key[i]) * FNV1A_PRIME; + } + return std::nullopt; + } + + // this mostly exists as a convenience + static std::optional 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(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` or `std::map` 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 table_rows_; + /*! tracks the row indices to make reverse lookups fast */ + std::shared_ptr 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(std::numeric_limits::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& keys) noexcept; + + /*! Lookup the value associated with the specified key + * + * This returns an empty optional if the key isn't known. + */ + std::optional find(const char* key) const noexcept; + std::optional find(std::string_view key) const noexcept; + /* + std::optional 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 std::string(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 + +#include // std::setfill, std::setw, std::hex +#include // 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 expected{{0, 0x811c9dc5ULL}}; + EXPECT_EQ(FNV1aHasher<>::calc(""), expected); + EXPECT_EQ(FNV1aHasher<>::calc(std::string_view("")), expected); +} + +TEST(tALLBiMapFNV1a, aString) +{ + std::optional expected{{1, 0xe40c292cULL}}; + EXPECT_EQ(FNV1aHasher<>::calc("a"), expected); + EXPECT_EQ(FNV1aHasher<>::calc(std::string_view("a")), expected); +} + +TEST(tALLBiMapFNV1a, foobarString) +{ + std::optional 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 expected{{6, 0xbf9cf968ULL}}; + EXPECT_EQ(FNV1aHasher::calc("foobar"), expected); + EXPECT_EQ(FNV1aHasher::calc(std::string_view("foobar")), expected); +} + +TEST(tALLBiMapFNV1a, TooLongString) +{ + constexpr int MaxKeyLen = 5; // <- shorter than the queried key + EXPECT_EQ(FNV1aHasher::calc("foobar"), std::nullopt); + EXPECT_EQ(FNV1aHasher::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 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{1}); + EXPECT_EQ(m.find(std::string("momentum_x")), std::optional{1}); + EXPECT_EQ(m.find(std::string_view("momentum_x")), std::optional{1}); + + // Equivalent Python/idiomatic C++: `6 == m["GasEnergy"]` + EXPECT_EQ(m.find("GasEnergy"), std::optional{6}); + EXPECT_EQ(m.find(std::string("GasEnergy")), std::optional{6}); + EXPECT_EQ(m.find(std::string_view("GasEnergy")), std::optional{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 keys = {"density", long_key}; + + ASSERT_DEATH({ utils::FrozenKeyIdxBiMap m(keys); }, "no key may have more than .* characters"); +} + +TEST(tALLBiMapDeathTest, EmptyKey) +{ + std::vector keys = {"density", ""}; + + ASSERT_DEATH({ utils::FrozenKeyIdxBiMap m(keys); }, "each key must hold at least 1 character"); +} + +TEST(tALLBiMapDeathTest, DuplicateKey) +{ + std::vector keys = {"density", "momentum_x", "density"}; + + ASSERT_DEATH({ utils::FrozenKeyIdxBiMap m(keys); }, "\"density\" key repeats"); +} \ No newline at end of file From 95ac2139eb3e219e236a3c79c2eae73bb63bba0f Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 21 Jan 2026 15:32:52 -0500 Subject: [PATCH 02/16] some more incremental progress --- src/grid/field_info.cpp | 61 +++++++++++++++++++++++++++++++++++++++++ src/grid/field_info.h | 11 ++++++++ src/grid/grid3D.cpp | 21 ++------------ src/grid/grid3D.h | 4 +++ 4 files changed, 79 insertions(+), 18 deletions(-) create mode 100644 src/grid/field_info.cpp create mode 100644 src/grid/field_info.h diff --git a/src/grid/field_info.cpp b/src/grid/field_info.cpp new file mode 100644 index 000000000..8ccff3154 --- /dev/null +++ b/src/grid/field_info.cpp @@ -0,0 +1,61 @@ +/*! \file + * Define machinery for accessing field information + */ + +#include "field_info.h" + +#include +#include + +#include "../utils/FrozenKeyIdxBiMap.h" +#include "grid_enum.h" + +/*! list of all field names + * + * This must remain synchronized with grid_enum.h + */ +static const char* field_names_[] = { + "density", "momentum_x", "momentum_y", "momentum_z", "Energy", + +#ifdef SCALAR + // Add scalars here, wrapped appropriately with ifdefs: + #ifdef BASIC_SCALAR + "basic_scalar", + #endif + + #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU) + "HI_density", "HII_density", "HeI_density", "HeII_density", "HeIII_density", "e_density", + #ifdef GRACKLE_METALS + "metal_density", + #endif + #endif + + #ifdef DUST + "dust_density", + #endif // DUST + +#endif // SCALAR + +#ifdef MHD + "magnetic_x", "magnetic_y", "magnetic_z", +#endif +#ifdef DE + "GasEnergy", +#endif +}; + +static constexpr std::size_t n_fields_ = sizeof(field_names_) / sizeof(const char*); + +static_assert(n_fields_ == grid_enum::num_fields, "field_names_ and grid_enum::num_fields are no longer synchronized"); + +utils::FrozenKeyIdxBiMap get_field_id_mapping() +{ + // convert n_fields to a vector of std::string + std::vector v; + v.reserve(n_fields_); + for (std::size_t i = 0; i < n_fields_; i++) { + v.push_back(std::string(field_names_[i])); + } + + return utils::FrozenKeyIdxBiMap(v); +} diff --git a/src/grid/field_info.h b/src/grid/field_info.h new file mode 100644 index 000000000..f5566dee4 --- /dev/null +++ b/src/grid/field_info.h @@ -0,0 +1,11 @@ +/*! \file + * Define machinery for accessing field information + */ + +#pragma once + +#include "../utils/FrozenKeyIdxBiMap.h" + +/*! Construct a bidirectional mapping between field name and field value + */ +utils::FrozenKeyIdxBiMap get_field_id_mapping(); \ No newline at end of file diff --git a/src/grid/grid3D.cpp b/src/grid/grid3D.cpp index 34b30a4ed..6b8ce1767 100644 --- a/src/grid/grid3D.cpp +++ b/src/grid/grid3D.cpp @@ -7,6 +7,7 @@ #include #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 @@ -112,24 +113,8 @@ 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 + field_name_map = get_field_id_mapping(); + H.n_fields = field_name_map.size(); int nx_in = P->nx; int ny_in = P->ny; diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h index 1cdaaab39..0cfd4b9f8 100644 --- a/src/grid/grid3D.h +++ b/src/grid/grid3D.h @@ -15,6 +15,7 @@ #include "../global/global.h" #include "../global/global_cuda.h" #include "../io/FnameTemplate.h" +#include "../utils/FrozenKeyIdxBiMap.h" #ifdef HDF5 #include @@ -297,6 +298,9 @@ class Grid3D * \brief Rotation struct for data projections */ struct Rotation R; + /*! Describes the mapping between field names and field indices */ + utils::FrozenKeyIdxBiMap field_name_map; + #ifdef GRAVITY // Object that contains data for gravity Grav3D Grav; From d0bb96b1c420ac38630ab61d7b1a32fe55413c51 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 22 Jan 2026 14:44:29 -0500 Subject: [PATCH 03/16] introduce FieldInfo --- src/grid/field_info.cpp | 71 +++++++++++++++++++++------ src/grid/field_info.h | 104 +++++++++++++++++++++++++++++++++++++++- src/grid/grid3D.cpp | 5 +- src/grid/grid3D.h | 4 +- 4 files changed, 163 insertions(+), 21 deletions(-) diff --git a/src/grid/field_info.cpp b/src/grid/field_info.cpp index 8ccff3154..6e710118b 100644 --- a/src/grid/field_info.cpp +++ b/src/grid/field_info.cpp @@ -10,52 +10,95 @@ #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; +}; + +} // anonymous namespace + /*! list of all field names * * This must remain synchronized with grid_enum.h */ -static const char* field_names_[] = { - "density", "momentum_x", "momentum_y", "momentum_z", "Energy", +static constexpr PropPack pack_arr_[] = {{"density", field::Kind::HYDRO}, {"momentum_x", field::Kind::HYDRO}, + {"momentum_y", field::Kind::HYDRO}, {"momentum_z", field::Kind::HYDRO}, + {"Energy", field::Kind::HYDRO}, #ifdef SCALAR // Add scalars here, wrapped appropriately with ifdefs: #ifdef BASIC_SCALAR - "basic_scalar", + {"basic_scalar", field::Kind::SCALAR}, #endif #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU) - "HI_density", "HII_density", "HeI_density", "HeII_density", "HeIII_density", "e_density", + {"HI_density", field::Kind::SCALAR}, {"HII_density", field::Kind::SCALAR}, + {"HeI_density", field::Kind::SCALAR}, {"HeII_density", field::Kind::SCALAR}, + {"HeIII_density", field::Kind::SCALAR}, {"e_density", field::Kind::SCALAR}, #ifdef GRACKLE_METALS - "metal_density", + {"metal_density", field::Kind::SCALAR}, #endif #endif #ifdef DUST - "dust_density", + {"dust_density", field::Kind::SCALAR}, #endif // DUST #endif // SCALAR #ifdef MHD - "magnetic_x", "magnetic_y", "magnetic_z", + {"magnetic_x", field::Kind::MAGNETIC}, {"magnetic_y", field::Kind::MAGNETIC}, + {"magnetic_z", field::Kind::MAGNETIC}, #endif #ifdef DE - "GasEnergy", + {"GasEnergy", field::Kind::HYDRO} #endif }; -static constexpr std::size_t n_fields_ = sizeof(field_names_) / sizeof(const char*); +static constexpr int n_fields_ = static_cast(sizeof(pack_arr_) / sizeof(PropPack)); -static_assert(n_fields_ == grid_enum::num_fields, "field_names_ and grid_enum::num_fields are no longer synchronized"); +static_assert(n_fields_ == grid_enum::num_fields, "pack_arr_ and grid_enum::num_fields are no longer synchronized"); -utils::FrozenKeyIdxBiMap get_field_id_mapping() +FieldInfo FieldInfo::create() { + FieldInfo out; + // convert n_fields to a vector of std::string std::vector v; v.reserve(n_fields_); for (std::size_t i = 0; i < n_fields_; i++) { - v.push_back(std::string(field_names_[i])); + v.push_back(std::string(pack_arr_[i].name)); + switch (pack_arr_[i].kind) { + case field::Kind::HYDRO: + out.hydro_field_ids_.push_back(i); + break; + case field::Kind::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"); + } } - - return utils::FrozenKeyIdxBiMap(v); + out.name_id_bimap_ = utils::FrozenKeyIdxBiMap(v); + return out; } + +const std::vector& FieldInfo::get_kind_ids_(field::Kind kind) const +{ + switch (kind) { + case field::Kind::HYDRO: + return hydro_field_ids_; + case field::Kind::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 index f5566dee4..d365921df 100644 --- a/src/grid/field_info.h +++ b/src/grid/field_info.h @@ -4,8 +4,108 @@ #pragma once +#include +#include +#include +#include + #include "../utils/FrozenKeyIdxBiMap.h" -/*! Construct a bidirectional mapping between field name and field value +namespace field +{ + +// note: HYDRO includes GasEnergy (if present) +enum class Kind { HYDRO, SCALAR, MAGNETIC }; + +/*! This is a "range" in the C++ 20 sense + * + * See \ref FieldInfo::get_id_range for an example */ -utils::FrozenKeyIdxBiMap get_field_id_mapping(); \ No newline at end of file +class IdRange +{ + const std::vector& id_vec_; + + public: + explicit IdRange(const std::vector& 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::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 hydro_field_ids_; + std::vector scalar_field_ids_; + std::vector magnetic_field_ids_; + + // 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& 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 field_id(const char* field_name) const { return name_id_bimap_.find(field_name); } + std::optional 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 field_name(int field_id) const + { + bool bad_id = (field_id < 0 || field_id >= n_fields()); + return bad_id ? std::nullopt : std::optional{name_id_bimap_.inverse_find(field_id)}; + } + + /*! Returns the number of fields */ + int n_fields() const { return static_cast(name_id_bimap_.size()); } + + /*! Returns the number of fields of a given category */ + int n_fields(field::Kind kind) const { return static_cast(get_kind_ids_(kind).size()); } + + /*! Returns the first field_id corresponding to a passive scalar (if there are any) */ + std::optional scalar_start() const + { + return scalar_field_ids_.empty() ? std::nullopt : std::optional{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 6b8ce1767..84d50a7d6 100644 --- a/src/grid/grid3D.cpp +++ b/src/grid/grid3D.cpp @@ -40,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; @@ -113,8 +113,7 @@ Real Grid3D::Calc_Inverse_Timestep() * \brief Initialize the grid. */ void Grid3D::Initialize(struct Parameters *P) { - field_name_map = get_field_id_mapping(); - H.n_fields = field_name_map.size(); + 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 0cfd4b9f8..deb21aa46 100644 --- a/src/grid/grid3D.h +++ b/src/grid/grid3D.h @@ -14,8 +14,8 @@ #include "../global/global.h" #include "../global/global_cuda.h" +#include "../grid/field_info.h" #include "../io/FnameTemplate.h" -#include "../utils/FrozenKeyIdxBiMap.h" #ifdef HDF5 #include @@ -299,7 +299,7 @@ class Grid3D struct Rotation R; /*! Describes the mapping between field names and field indices */ - utils::FrozenKeyIdxBiMap field_name_map; + FieldInfo field_info; #ifdef GRAVITY // Object that contains data for gravity From 4c0c701caf331b8efdca7f1facbf070e4c85b7bc Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 22 Jan 2026 14:51:50 -0500 Subject: [PATCH 04/16] eliminate unused output_ionization variable --- src/io/io.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index 93178fc09..597fc6689 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -1322,7 +1322,7 @@ void Grid3D::Write_Grid_HDF5(hid_t file_id) #endif // OUTPUT_MOMENTUM #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU) - bool output_metals, output_electrons, output_full_ionization; + bool output_metals, output_electrons; #ifdef OUTPUT_METALS output_metals = true; #else // not OUTPUT_METALS @@ -1333,11 +1333,6 @@ void Grid3D::Write_Grid_HDF5(hid_t file_id) #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 From d1952cb533033b843859dc4029eaefb6840c42e9 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 27 Jan 2026 09:17:34 -0500 Subject: [PATCH 05/16] a bunch more progress --- src/grid/field_info.cpp | 35 ++++++---- src/grid/field_info.h | 14 ++++ src/grid/grid3D.h | 8 ++- src/io/FieldWriter.cpp | 95 +++++++++++++++++++++++++++ src/io/FieldWriter.h | 47 +++++++++++++ src/io/WriterManager.cpp | 6 +- src/io/WriterManager.h | 3 +- src/io/io.cpp | 113 ++++++-------------------------- src/io/io.h | 4 +- src/main.cpp | 2 +- src/utils/FrozenKeyIdxBiMap.cpp | 12 ++-- src/utils/FrozenKeyIdxBiMap.h | 2 +- 12 files changed, 223 insertions(+), 118 deletions(-) create mode 100644 src/io/FieldWriter.cpp create mode 100644 src/io/FieldWriter.h diff --git a/src/grid/field_info.cpp b/src/grid/field_info.cpp index 6e710118b..4d15d7356 100644 --- a/src/grid/field_info.cpp +++ b/src/grid/field_info.cpp @@ -16,6 +16,7 @@ namespace struct PropPack { const char* name; field::Kind kind; + field::IOBuf io_buf; }; } // anonymous namespace @@ -24,37 +25,44 @@ struct PropPack { * * This must remain synchronized with grid_enum.h */ -static constexpr PropPack pack_arr_[] = {{"density", field::Kind::HYDRO}, {"momentum_x", field::Kind::HYDRO}, - {"momentum_y", field::Kind::HYDRO}, {"momentum_z", field::Kind::HYDRO}, - {"Energy", field::Kind::HYDRO}, +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 - // Add scalars here, wrapped appropriately with ifdefs: #ifdef BASIC_SCALAR - {"basic_scalar", field::Kind::SCALAR}, + // we use the name "scalar0" for better consistency with the name recorded during IO + {"scalar0", field::Kind::SCALAR, field::IOBuf::DEVICE}, #endif #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU) - {"HI_density", field::Kind::SCALAR}, {"HII_density", field::Kind::SCALAR}, - {"HeI_density", field::Kind::SCALAR}, {"HeII_density", field::Kind::SCALAR}, - {"HeIII_density", field::Kind::SCALAR}, {"e_density", field::Kind::SCALAR}, + {"HI_density", field::Kind::SCALAR, field::IOBuf::HOST}, + {"HII_density", field::Kind::SCALAR, field::IOBuf::HOST}, + {"HeI_density", field::Kind::SCALAR, field::IOBuf::HOST}, + {"HeII_density", field::Kind::SCALAR, field::IOBuf::HOST}, + {"HeIII_density", field::Kind::SCALAR, field::IOBuf::HOST}, + {"e_density", field::Kind::SCALAR, field::IOBuf::HOST}, #ifdef GRACKLE_METALS - {"metal_density", field::Kind::SCALAR}, + {"metal_density", field::Kind::SCALAR, field::IOBuf::HOST}, #endif #endif #ifdef DUST - {"dust_density", field::Kind::SCALAR}, + {"dust_density", field::Kind::SCALAR, field::IOBuf::DEVICE}, #endif // DUST #endif // SCALAR #ifdef MHD - {"magnetic_x", field::Kind::MAGNETIC}, {"magnetic_y", field::Kind::MAGNETIC}, - {"magnetic_z", field::Kind::MAGNETIC}, + {"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} + {"GasEnergy", field::Kind::HYDRO, eld::IOBuf::DEVICE} #endif }; @@ -71,6 +79,7 @@ FieldInfo FieldInfo::create() 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); diff --git a/src/grid/field_info.h b/src/grid/field_info.h index d365921df..c8f4d167e 100644 --- a/src/grid/field_info.h +++ b/src/grid/field_info.h @@ -17,6 +17,9 @@ namespace field // note: HYDRO includes GasEnergy (if present) enum class Kind { HYDRO, SCALAR, MAGNETIC }; +/*! Specifies which buffer to use for IO */ +enum class IOBuf { HOST, DEVICE }; + /*! This is a "range" in the C++ 20 sense * * See \ref FieldInfo::get_id_range for an example @@ -46,6 +49,7 @@ class FieldInfo std::vector hydro_field_ids_; std::vector scalar_field_ids_; std::vector magnetic_field_ids_; + std::vector io_buf_; // We make the default-constructor private to force the use of the factory method FieldInfo() = default; @@ -86,6 +90,16 @@ class FieldInfo return bad_id ? std::nullopt : std::optional{name_id_bimap_.inverse_find(field_id)}; } + /*! 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 io_buf(int field_id) const + { + bool bad_id = (field_id < 0 || field_id >= n_fields()); + return bad_id ? std::nullopt : std::optional{io_buf_[field_id]}; + } + /*! Returns the number of fields */ int n_fields() const { return static_cast(name_id_bimap_.size()); } diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h index deb21aa46..d88b4c119 100644 --- a/src/grid/grid3D.h +++ b/src/grid/grid3D.h @@ -51,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*/ @@ -503,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_arr, int n_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/io/FieldWriter.cpp b/src/io/FieldWriter.cpp new file mode 100644 index 000000000..d0e7cb8d4 --- /dev/null +++ b/src/io/FieldWriter.cpp @@ -0,0 +1,95 @@ +/*! + * \file + * Implements the FieldWriter type + */ + +#include "../io/FieldWriter.h" + +#include +#include +#include +#include + +#include "../grid/field_info.h" +#include "../io/io.h" +#include "../utils/error_handling.h" + +namespace io +{ + +// 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) +{ + std::vector& vec = this->h5_dataset_spec_; + auto add_dataset_entry = [&vec, &field_info](const char* name, WriteCond cond) { + std::optional 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}); + }; + + add_dataset_entry("density", WriteCond::ALWAYS); + add_dataset_entry("momentum_x", MOMENTUM_CONDITION); + add_dataset_entry("momentum_y", MOMENTUM_CONDITION); + add_dataset_entry("momentum_z", MOMENTUM_CONDITION); + add_dataset_entry("Energy", ENERGY_CONDITION); +#ifdef DE + add_dataset_entry("GasEnergy", ENERGY_CONDITION); +#endif + + for (int field_id : field_info.get_id_range(field::Kind::SCALAR)) { + std::string name = field_info.field_name(field_id).value(); + if (name == "e_density") { + add_dataset_entry(name.c_str(), ELECTRONS_CONDITION); + } else if (name == "metal_density") { + add_dataset_entry(name.c_str(), METALS_CONDITION); + } else { + add_dataset_entry(name.c_str(), WriteCond::ALWAYS); + } + } + + // For now, I'm intentionally ignoring the remaining assorted outputs (e.g. magnetic + // fields, temperature, gravitational potential). That stuff is still handled very + // manually) +} + +void FieldWriter::operator()(Grid3D& G, Parameters P, int nfile, const FnameTemplate& fname_template) const +{ + Output_Data(G, P, nfile, fname_template, h5_dataset_spec_); +} + +} // 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..0a69a6eb5 --- /dev/null +++ b/src/io/FieldWriter.h @@ -0,0 +1,47 @@ +/*! + * \file + * Declares the FieldWriter type + */ + +#pragma once + +#include +#include + +#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 + +namespace io +{ + +enum struct WriteCond { ALWAYS, REQUIRE_COMPLETE_DATA }; + +struct DatasetSpec { + 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; +}; + +/*! \brief A callable that writes general grid data + */ +class FieldWriter +{ + std::vector 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; +}; + +} // namespace io \ No newline at end of file diff --git a/src/io/WriterManager.cpp b/src/io/WriterManager.cpp index 02631b317..8f62f6745 100644 --- a/src/io/WriterManager.cpp +++ b/src/io/WriterManager.cpp @@ -10,10 +10,12 @@ #include #include "../gravity/grav3D.h" +#include "../io/FieldWriter.h" // FieldWriter #include "../io/ParameterMap.h" // define ParameterMap #include "../io/io.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) { // in the future, the goal is to read directly from ParameterMap (so we can stop storing // some of the relevant variables in Parameters) @@ -21,7 +23,7 @@ io::WriterManager::WriterManager(const Parameters& P, ParameterMap& pmap) : fnam #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) 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 #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 597fc6689..3ec6439eb 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -145,7 +145,8 @@ void Write_Data(Grid3D &G, struct Parameters P, int nfile, const io::WriterManag } /* Output the grid data to file. */ -void Output_Data(Grid3D &G, struct Parameters P, int nfile, const FnameTemplate &fname_template) +void Output_Data(Grid3D &G, struct Parameters P, int nfile, const FnameTemplate &fname_template, + const std::vector &h5_dataset_spec) { // create the filename std::string filename = fname_template.format_fname(nfile, ""); @@ -184,7 +185,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.data(), h5_dataset_spec.size()); // close the file status = H5Fclose(file_id); @@ -1298,7 +1299,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_arr, int n_dataset_spec) { int i, j, k, id, buf_id; hid_t dataset_id, dataspace_id; @@ -1306,36 +1307,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; - #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 - - #endif // COOLING_GRACKLE or CHEMISTRY_GPU - // Allocate necessary buffers int nx_dset = H.nx_real; int ny_dset = H.ny_real; @@ -1349,72 +1320,26 @@ 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 (int i = 0; i < n_dataset_spec; i++) { + const io::DatasetSpec cur_spec = h5_dataset_spec_arr[i]; + 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) { diff --git a/src/io/io.h b/src/io/io.h index 32e03ff8b..4c39787a2 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" @@ -27,7 +28,8 @@ void Print_Stats(Grid3D& G); 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); +void Output_Data(Grid3D& G, struct Parameters P, int nfile, const FnameTemplate& fname_template, + const std::vector& h5_dataset_spec); /* Output the grid data to file as 32-bit floats. */ void Output_Float32(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 index ddfc709dc..61131dc52 100644 --- a/src/utils/FrozenKeyIdxBiMap.cpp +++ b/src/utils/FrozenKeyIdxBiMap.cpp @@ -2,7 +2,9 @@ * Implementation for FrozenKeyIdxBiMap */ -#include "FrozenKeyIdxBiMap.h" +#include "../utils/FrozenKeyIdxBiMap.h" + +#include "../utils/error_handling.h" /*! Defines a bunch of extra details for implementing \ref FrozenKeyIdxBiMap */ namespace utils @@ -103,7 +105,8 @@ FrozenKeyIdxBiMap::FrozenKeyIdxBiMap(const std::vector& keys) noexc } // allocate hash table (each row is default constructed) - bimap_detail::Row* rows = new bimap_detail::Row[capacity]; + 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++) { @@ -114,8 +117,9 @@ FrozenKeyIdxBiMap::FrozenKeyIdxBiMap(const std::vector& keys) noexc delete[] rows; }; this->table_rows_ = std::shared_ptr(rows, table_deleter); - uint16_t* ordered_row_indices = new uint16_t[n_keys]; - this->ordered_row_indices_ = std::shared_ptr(ordered_row_indices, std::default_delete()); + 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(ordered_row_indices, std::default_delete()); this->length_ = n_keys; this->capacity_ = capacity; diff --git a/src/utils/FrozenKeyIdxBiMap.h b/src/utils/FrozenKeyIdxBiMap.h index ac109d504..d7d660aa9 100644 --- a/src/utils/FrozenKeyIdxBiMap.h +++ b/src/utils/FrozenKeyIdxBiMap.h @@ -179,7 +179,7 @@ class FrozenKeyIdxBiMap std::string inverse_find(int index) const { uint16_t row_index = ordered_row_indices_.get()[index]; - return std::string(table_rows_.get()[row_index].key); + return {table_rows_.get()[row_index].key}; } /*! return the number of keys in the map */ From 09072141e0379c3cba795c2c6b799674d56b4c32 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 28 Jan 2026 17:21:52 -0500 Subject: [PATCH 06/16] remove a few unused commented out lines --- src/io/io.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index 3ec6439eb..11a2405a6 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -30,9 +30,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); From 98e89c5cae5d7e18342476375d39968523ed3780 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 28 Jan 2026 17:26:21 -0500 Subject: [PATCH 07/16] minor typo fix --- src/grid/field_info.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/grid/field_info.cpp b/src/grid/field_info.cpp index 4d15d7356..0ccb7749b 100644 --- a/src/grid/field_info.cpp +++ b/src/grid/field_info.cpp @@ -62,7 +62,7 @@ static constexpr PropPack pack_arr_[] = { {"magnetic_z", field::Kind::MAGNETIC, field::IOBuf::DEVICE}, #endif #ifdef DE - {"GasEnergy", field::Kind::HYDRO, eld::IOBuf::DEVICE} + {"GasEnergy", field::Kind::HYDRO, field::IOBuf::DEVICE} #endif }; From 2f526b9a97b5940e85f07f99200effa56bb5e3d1 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 28 Jan 2026 17:30:01 -0500 Subject: [PATCH 08/16] Add a clarifying comment --- src/grid/grid_enum.h | 3 +++ 1 file changed, 3 insertions(+) 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 { From abe8892d604502d979677de786a3b8c8c92f64ff Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 3 Feb 2026 09:08:56 -0500 Subject: [PATCH 09/16] rename field::Kind::{SCALAR->PASSIVE_SCALAR} --- src/grid/field_info.cpp | 22 +++++++++++----------- src/grid/field_info.h | 27 ++++++++++++++++++++++++++- src/io/FieldWriter.cpp | 2 +- 3 files changed, 38 insertions(+), 13 deletions(-) diff --git a/src/grid/field_info.cpp b/src/grid/field_info.cpp index 0ccb7749b..b79dc2119 100644 --- a/src/grid/field_info.cpp +++ b/src/grid/field_info.cpp @@ -35,23 +35,23 @@ static constexpr PropPack pack_arr_[] = { #ifdef SCALAR #ifdef BASIC_SCALAR // we use the name "scalar0" for better consistency with the name recorded during IO - {"scalar0", field::Kind::SCALAR, field::IOBuf::DEVICE}, + {"scalar0", field::Kind::PASSIVE_SCALAR, field::IOBuf::DEVICE}, #endif #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU) - {"HI_density", field::Kind::SCALAR, field::IOBuf::HOST}, - {"HII_density", field::Kind::SCALAR, field::IOBuf::HOST}, - {"HeI_density", field::Kind::SCALAR, field::IOBuf::HOST}, - {"HeII_density", field::Kind::SCALAR, field::IOBuf::HOST}, - {"HeIII_density", field::Kind::SCALAR, field::IOBuf::HOST}, - {"e_density", field::Kind::SCALAR, field::IOBuf::HOST}, + {"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::SCALAR, field::IOBuf::HOST}, + {"metal_density", field::Kind::PASSIVE_SCALAR, field::IOBuf::HOST}, #endif #endif #ifdef DUST - {"dust_density", field::Kind::SCALAR, field::IOBuf::DEVICE}, + {"dust_density", field::Kind::PASSIVE_SCALAR, field::IOBuf::DEVICE}, #endif // DUST #endif // SCALAR @@ -84,7 +84,7 @@ FieldInfo FieldInfo::create() case field::Kind::HYDRO: out.hydro_field_ids_.push_back(i); break; - case field::Kind::SCALAR: + case field::Kind::PASSIVE_SCALAR: out.scalar_field_ids_.push_back(i); break; case field::Kind::MAGNETIC: @@ -103,7 +103,7 @@ const std::vector& FieldInfo::get_kind_ids_(field::Kind kind) const switch (kind) { case field::Kind::HYDRO: return hydro_field_ids_; - case field::Kind::SCALAR: + case field::Kind::PASSIVE_SCALAR: return scalar_field_ids_; case field::Kind::MAGNETIC: return magnetic_field_ids_; diff --git a/src/grid/field_info.h b/src/grid/field_info.h index c8f4d167e..8af559874 100644 --- a/src/grid/field_info.h +++ b/src/grid/field_info.h @@ -15,11 +15,13 @@ namespace field { // note: HYDRO includes GasEnergy (if present) -enum class Kind { HYDRO, SCALAR, MAGNETIC }; +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 @@ -90,6 +92,29 @@ class FieldInfo return bad_id ? std::nullopt : std::optional{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. 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 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{false}; + } + } + return std::optional{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. diff --git a/src/io/FieldWriter.cpp b/src/io/FieldWriter.cpp index d0e7cb8d4..56e254724 100644 --- a/src/io/FieldWriter.cpp +++ b/src/io/FieldWriter.cpp @@ -71,7 +71,7 @@ FieldWriter::FieldWriter(ParameterMap& pmap, const FieldInfo& field_info) add_dataset_entry("GasEnergy", ENERGY_CONDITION); #endif - for (int field_id : field_info.get_id_range(field::Kind::SCALAR)) { + 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") { add_dataset_entry(name.c_str(), ELECTRONS_CONDITION); From d7aa63ea31da81f4d400cafad3ccc3c64b5e1fb7 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 3 Feb 2026 09:15:38 -0500 Subject: [PATCH 10/16] convert Output_Data into the callable method of FieldWriter --- src/io/FieldWriter.cpp | 5 ----- src/io/io.cpp | 7 +++---- src/io/io.h | 4 ---- 3 files changed, 3 insertions(+), 13 deletions(-) diff --git a/src/io/FieldWriter.cpp b/src/io/FieldWriter.cpp index 56e254724..69cab8e43 100644 --- a/src/io/FieldWriter.cpp +++ b/src/io/FieldWriter.cpp @@ -87,9 +87,4 @@ FieldWriter::FieldWriter(ParameterMap& pmap, const FieldInfo& field_info) // manually) } -void FieldWriter::operator()(Grid3D& G, Parameters P, int nfile, const FnameTemplate& fname_template) const -{ - Output_Data(G, P, nfile, fname_template, h5_dataset_spec_); -} - } // namespace io \ No newline at end of file diff --git a/src/io/io.cpp b/src/io/io.cpp index 11a2405a6..9e29054e2 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -14,6 +14,7 @@ #include #endif // HDF5 #include "../grid/grid3D.h" +#include "../io/FieldWriter.h" #include "../io/WriterManager.h" #include "../io/io.h" #include "../utils/cuda_utilities.h" @@ -141,9 +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, - const std::vector &h5_dataset_spec) +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, ""); @@ -182,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, h5_dataset_spec.data(), h5_dataset_spec.size()); + G.Write_Grid_HDF5(file_id, h5_dataset_spec_.data(), h5_dataset_spec_.size()); // close the file status = H5Fclose(file_id); diff --git a/src/io/io.h b/src/io/io.h index 4c39787a2..68c63f533 100644 --- a/src/io/io.h +++ b/src/io/io.h @@ -27,10 +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, - const std::vector& h5_dataset_spec); - /* Output the grid data to file as 32-bit floats. */ void Output_Float32(Grid3D& G, struct Parameters P, int nfile, const FnameTemplate& fname_template); From f4fedb5830b669bcc93eb3dc84a3fab3c9f4c238 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 3 Feb 2026 09:30:43 -0500 Subject: [PATCH 11/16] light refactoring --- src/grid/grid3D.h | 2 +- src/io/FieldWriter.cpp | 15 ++++++++++----- src/io/FieldWriter.h | 20 ++++++++++++++++---- src/io/io.cpp | 7 +++---- 4 files changed, 30 insertions(+), 14 deletions(-) diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h index d88b4c119..b0028edf0 100644 --- a/src/grid/grid3D.h +++ b/src/grid/grid3D.h @@ -509,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, const io::DatasetSpec *h5_dataset_spec_arr, int n_dataset_spec); + 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/io/FieldWriter.cpp b/src/io/FieldWriter.cpp index 69cab8e43..75263ed89 100644 --- a/src/io/FieldWriter.cpp +++ b/src/io/FieldWriter.cpp @@ -52,8 +52,8 @@ static constexpr WriteCond ELECTRONS_CONDITION = WriteCond::REQUIRE_COMPLETE_DAT FieldWriter::FieldWriter(ParameterMap& pmap, const FieldInfo& field_info) { - std::vector& vec = this->h5_dataset_spec_; - auto add_dataset_entry = [&vec, &field_info](const char* name, WriteCond cond) { + std::vector& vec = this->h5_dataset_spec_.cc_dataset_entries; + auto add_dataset_entry = [&vec, &field_info](const char* name, WriteCond cond) { std::optional 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); @@ -82,9 +82,14 @@ FieldWriter::FieldWriter(ParameterMap& pmap, const FieldInfo& field_info) } } - // For now, I'm intentionally ignoring the remaining assorted outputs (e.g. magnetic - // fields, temperature, gravitational potential). That stuff is still handled very - // manually) +#ifdef MHD + h5_dataset_spec_.mhd_condition = std::optional{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) } } // namespace io \ No newline at end of file diff --git a/src/io/FieldWriter.h b/src/io/FieldWriter.h index 0a69a6eb5..fc993aade 100644 --- a/src/io/FieldWriter.h +++ b/src/io/FieldWriter.h @@ -5,21 +5,23 @@ #pragma once +#include #include #include #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 "../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 DatasetSpec { +struct DatasetSpecEntry { int field_id; /// the dataset name. By convention, this is prefixed with a "/" std::string name; @@ -29,11 +31,21 @@ struct DatasetSpec { WriteCond condition; }; +struct DatasetSpec { + /// describes properties about dataset creation for ordinary cell-centered fields + std::vector 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 mhd_condition; +}; + /*! \brief A callable that writes general grid data */ class FieldWriter { - std::vector h5_dataset_spec_; + DatasetSpec h5_dataset_spec_; public: FieldWriter() = delete; diff --git a/src/io/io.cpp b/src/io/io.cpp index 9e29054e2..1c51afc6f 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -181,7 +181,7 @@ void io::FieldWriter::operator()(Grid3D &G, Parameters P, int nfile, const Fname G.Write_Header_HDF5(file_id); // write the conserved variables to the output file - G.Write_Grid_HDF5(file_id, h5_dataset_spec_.data(), h5_dataset_spec_.size()); + G.Write_Grid_HDF5(file_id, h5_dataset_spec_); // close the file status = H5Fclose(file_id); @@ -1295,7 +1295,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, const io::DatasetSpec *h5_dataset_spec_arr, int n_dataset_spec) +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; @@ -1316,8 +1316,7 @@ void Grid3D::Write_Grid_HDF5(hid_t file_id, const io::DatasetSpec *h5_dataset_sp dataset_buffer = (Real *)malloc(buffer_size * sizeof(Real)); // Start writing fields - for (int i = 0; i < n_dataset_spec; i++) { - const io::DatasetSpec cur_spec = h5_dataset_spec_arr[i]; + 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; } From ebad964a7e921862ab28e30f53ea21866db6ed82 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 3 Feb 2026 10:33:24 -0500 Subject: [PATCH 12/16] slightly adjust handling of Magnetic fields --- src/io/io.cpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index 1c51afc6f..2c8886a5c 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -1345,16 +1345,17 @@ void Grid3D::Write_Grid_HDF5(hid_t file_id, const io::DatasetSpec &h5_dataset_sp 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); From 47488aeb97f37fdc706741ac4935dbaca667d1cd Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 2 Feb 2026 12:52:33 -0500 Subject: [PATCH 13/16] an initial commit --- src/io/FieldWriter.cpp | 4 ++++ src/io/FieldWriter.h | 21 +++++++++++++++++++++ src/io/WriterManager.cpp | 2 +- src/io/io.cpp | 2 +- src/io/io.h | 3 --- 5 files changed, 27 insertions(+), 5 deletions(-) diff --git a/src/io/FieldWriter.cpp b/src/io/FieldWriter.cpp index 75263ed89..475e32113 100644 --- a/src/io/FieldWriter.cpp +++ b/src/io/FieldWriter.cpp @@ -92,4 +92,8 @@ FieldWriter::FieldWriter(ParameterMap& pmap, const FieldInfo& field_info) // temperature, gravitational potential). That stuff is still handled very manually) } +F32FieldWriter::F32FieldWriter(ParameterMap& pmap, const FieldInfo &field_info) +{ +} + } // namespace io \ No newline at end of file diff --git a/src/io/FieldWriter.h b/src/io/FieldWriter.h index fc993aade..ad320abc7 100644 --- a/src/io/FieldWriter.h +++ b/src/io/FieldWriter.h @@ -5,6 +5,7 @@ #pragma once +#include #include #include #include @@ -42,6 +43,8 @@ struct DatasetSpec { }; /*! \brief A callable that writes general grid data + * + * \todo Maybe work to consolidate this with F32FieldWriter */ class FieldWriter { @@ -56,4 +59,22 @@ class FieldWriter 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::array write_mag_xyz = {false, false, false}; + std::vector main_fields; + + 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 8f62f6745..fcba3a1d4 100644 --- a/src/io/WriterManager.cpp +++ b/src/io/WriterManager.cpp @@ -29,7 +29,7 @@ io::WriterManager::WriterManager(const Parameters& P, ParameterMap& pmap, const // 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}); + packs_.push_back(io::detail::WriterPack{"hydro-f32", n_hydro, {io::F32FieldWriter(pmap, field_info)}}); } #endif diff --git a/src/io/io.cpp b/src/io/io.cpp index 2c8886a5c..329434c8c 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -211,7 +211,7 @@ void io::FieldWriter::operator()(Grid3D &G, Parameters P, int nfile, const Fname #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; diff --git a/src/io/io.h b/src/io/io.h index 68c63f533..179312ccd 100644 --- a/src/io/io.h +++ b/src/io/io.h @@ -27,9 +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 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); From ada83bac56a1561179bde22b582542560dd9689b Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 2 Feb 2026 13:42:05 -0500 Subject: [PATCH 14/16] shift extra OUT_FLOAT_32 checks out of callback --- src/io/WriterManager.cpp | 23 ++++++++++++++++++++--- src/io/io.cpp | 11 ----------- 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/src/io/WriterManager.cpp b/src/io/WriterManager.cpp index fcba3a1d4..2992a0b49 100644 --- a/src/io/WriterManager.cpp +++ b/src/io/WriterManager.cpp @@ -6,6 +6,8 @@ #include "../io/WriterManager.h" #include +#include +#include // std::lcm #include #include @@ -13,23 +15,38 @@ #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, 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, {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, {io::F32FieldWriter(pmap, field_info)}}); + // 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::max()}, + "the lcm of n_hydro and n_out_float32 can't be represented by an int"); + int cadence = static_cast(lcm); + + packs_.push_back(io::detail::WriterPack{"hydro-f32", cadence, {io::F32FieldWriter(pmap, field_info)}}); } #endif diff --git a/src/io/io.cpp b/src/io/io.cpp index 329434c8c..cac983684 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -215,17 +215,6 @@ void io::F32FieldWriter::operator()(Grid3D &G, Parameters P, int nfile, const Fn { #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"); From 22cac7a8a438d401fbb26c12d8b4c54a6d8c2d39 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 2 Feb 2026 14:21:31 -0500 Subject: [PATCH 15/16] Factor lambda function out of FieldWriter::FieldWriter --- src/io/FieldWriter.cpp | 56 +++++++++++++++++++++++++----------------- 1 file changed, 34 insertions(+), 22 deletions(-) diff --git a/src/io/FieldWriter.cpp b/src/io/FieldWriter.cpp index 475e32113..dfaf0c567 100644 --- a/src/io/FieldWriter.cpp +++ b/src/io/FieldWriter.cpp @@ -17,6 +17,27 @@ 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& vec; + const FieldInfo& field_info; + + void add_entry(const char* name, WriteCond cond) + { + std::optional 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 @@ -52,33 +73,26 @@ static constexpr WriteCond ELECTRONS_CONDITION = WriteCond::REQUIRE_COMPLETE_DAT FieldWriter::FieldWriter(ParameterMap& pmap, const FieldInfo& field_info) { - std::vector& vec = this->h5_dataset_spec_.cc_dataset_entries; - auto add_dataset_entry = [&vec, &field_info](const char* name, WriteCond cond) { - std::optional 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}); - }; - - add_dataset_entry("density", WriteCond::ALWAYS); - add_dataset_entry("momentum_x", MOMENTUM_CONDITION); - add_dataset_entry("momentum_y", MOMENTUM_CONDITION); - add_dataset_entry("momentum_z", MOMENTUM_CONDITION); - add_dataset_entry("Energy", ENERGY_CONDITION); + // construct DsetSpecListBuilder_ to append entries to `this->h5_dataset_spec_` + 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 - add_dataset_entry("GasEnergy", ENERGY_CONDITION); + 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") { - add_dataset_entry(name.c_str(), ELECTRONS_CONDITION); + dsentry_l_builder.add_entry(name.c_str(), ELECTRONS_CONDITION); } else if (name == "metal_density") { - add_dataset_entry(name.c_str(), METALS_CONDITION); + dsentry_l_builder.add_entry(name.c_str(), METALS_CONDITION); } else { - add_dataset_entry(name.c_str(), WriteCond::ALWAYS); + dsentry_l_builder.add_entry(name.c_str(), WriteCond::ALWAYS); } } @@ -92,8 +106,6 @@ FieldWriter::FieldWriter(ParameterMap& pmap, const FieldInfo& field_info) // temperature, gravitational potential). That stuff is still handled very manually) } -F32FieldWriter::F32FieldWriter(ParameterMap& pmap, const FieldInfo &field_info) -{ -} +F32FieldWriter::F32FieldWriter(ParameterMap& pmap, const FieldInfo& field_info) {} } // namespace io \ No newline at end of file From 346444a4c7d1127514bdd815a3182ad3281b2888 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 3 Feb 2026 10:58:32 -0500 Subject: [PATCH 16/16] incremental commit --- src/global/global.cpp | 14 -------------- src/global/global.h | 26 +++++++++----------------- src/io/FieldWriter.cpp | 17 +++++++++++++++-- src/io/FieldWriter.h | 4 +--- src/io/io.cpp | 35 +++++++++++------------------------ 5 files changed, 36 insertions(+), 60 deletions(-) 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/io/FieldWriter.cpp b/src/io/FieldWriter.cpp index dfaf0c567..6fb45551c 100644 --- a/src/io/FieldWriter.cpp +++ b/src/io/FieldWriter.cpp @@ -73,7 +73,7 @@ static constexpr WriteCond ELECTRONS_CONDITION = WriteCond::REQUIRE_COMPLETE_DAT FieldWriter::FieldWriter(ParameterMap& pmap, const FieldInfo& field_info) { - // construct DsetSpecListBuilder_ to append entries to `this->h5_dataset_spec_` + // 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); @@ -106,6 +106,19 @@ FieldWriter::FieldWriter(ParameterMap& pmap, const FieldInfo& field_info) // temperature, gravitational potential). That stuff is still handled very manually) } -F32FieldWriter::F32FieldWriter(ParameterMap& pmap, const FieldInfo& field_info) {} +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 index ad320abc7..30427175b 100644 --- a/src/io/FieldWriter.h +++ b/src/io/FieldWriter.h @@ -5,7 +5,6 @@ #pragma once -#include #include #include #include @@ -65,8 +64,7 @@ class FieldWriter */ class F32FieldWriter { - std::array write_mag_xyz = {false, false, false}; - std::vector main_fields; + std::vector cc_dataset_entries; public: F32FieldWriter() = delete; diff --git a/src/io/io.cpp b/src/io/io.cpp index cac983684..481748da3 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -252,32 +252,19 @@ void io::F32FieldWriter::operator()(Grid3D &G, Parameters P, int nfile, const Fn cuda_utilities::DeviceVector static device_dataset_vector{buffer_size}; auto *dataset_buffer = (float *)malloc(buffer_size * sizeof(float)); - if (P.out_float32_density > 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_density, "/density"); + device_dataset_vector.data(), ptr, cur_spec.name.c_str()); } - 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) { - 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"); - } - #endif // DE + + // TODO: finish transitioning logic #ifdef MHD // TODO (by Alwin, for anyone) : Repair output format if needed and remove these chprintfs when appropriate