hdf_common.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2022-2023 Dennis Gläser <dennis.glaeser@iws.uni-stuttgart.de>
// SPDX-License-Identifier: MIT

#ifndef GRIDFORMAT_VTK_HDF_COMMON_HPP_
#define GRIDFORMAT_VTK_HDF_COMMON_HPP_
 
#include <vector>
#include <cstddef>
#include <numeric>
#include <string>
#include <optional>
#include <ranges>
#include <array>
#include <exception>
 
#include <gridformat/common/hdf5.hpp>
#include <gridformat/common/field.hpp>
#include <gridformat/common/precision.hpp>
#include <gridformat/common/lazy_field.hpp>
#include <gridformat/common/range_field.hpp>
#include <gridformat/common/field_transformations.hpp>
#include <gridformat/common/multi_string.hpp>
#include <gridformat/common/string_conversion.hpp>
#include <gridformat/common/exceptions.hpp>
#include <gridformat/parallel/communication.hpp>
 
namespace GridFormat {
 
namespace VTK {

struct HDFTransientOptions {
    bool static_grid = false; 
    bool static_meta_data = true; 
};
 
}  // namespace VTK
 
 
namespace VTKHDF {

struct IOContext {
    const int my_rank;
    const int num_ranks;
    const bool is_parallel;
    const std::vector<std::size_t> rank_cells;
    const std::vector<std::size_t> rank_points;
    const std::size_t num_cells_total;
    const std::size_t num_points_total;
    const std::size_t my_cell_offset;
    const std::size_t my_point_offset;
 
    IOContext(int _my_rank,
              int _num_ranks,
              std::vector<std::size_t> _rank_cells,
              std::vector<std::size_t> _rank_points)
    : my_rank{_my_rank}
    , num_ranks{_num_ranks}
    , is_parallel{_num_ranks > 1}
    , rank_cells{std::move(_rank_cells)}
    , rank_points{std::move(_rank_points)}
    , num_cells_total{_accumulate(rank_cells)}
    , num_points_total{_accumulate(rank_points)}
    , my_cell_offset{_accumulate_rank_offset(rank_cells)}
    , my_point_offset{_accumulate_rank_offset(rank_points)} {
        if (my_rank >= num_ranks)
            throw ValueError(as_error("Given rank is not within communicator size"));
        if (num_ranks != static_cast<int>(rank_cells.size()))
            throw ValueError(as_error("Cells vector does not match communicator size"));
        if (num_ranks != static_cast<int>(rank_points.size()))
            throw ValueError(as_error("Points vector does not match communicator size"));
    }
 
    template<Concepts::Grid Grid, Concepts::Communicator Communicator>
    static IOContext from(const Grid& grid, const Communicator& comm, int root_rank = 0) {
        const int size = Parallel::size(comm);
        const int rank = Parallel::rank(comm);
        const std::size_t num_points = number_of_points(grid);
        const std::size_t num_cells = number_of_cells(grid);
        if (size == 1)
            return IOContext{rank, size, std::vector{num_cells}, std::vector{num_points}};
 
        const auto all_num_points = Parallel::gather(comm, num_points, root_rank);
        const auto all_num_cells = Parallel::gather(comm, num_cells, root_rank);
        const auto my_all_num_points = Parallel::broadcast(comm, all_num_points, root_rank);
        const auto my_all_num_cells = Parallel::broadcast(comm, all_num_cells, root_rank);
        return IOContext{rank, size, my_all_num_cells, my_all_num_points};
    }
 
 private:
    std::size_t _accumulate(const std::vector<std::size_t>& in) const {
        return std::accumulate(in.begin(), in.end(), std::size_t{0});
    }
 
    std::size_t _accumulate_rank_offset(const std::vector<std::size_t>& in) const {
        if (in.size() <= static_cast<std::size_t>(my_rank))
            throw ValueError("Rank-vector length must be equal to number of ranks");
        return std::accumulate(in.begin(), std::next(in.begin(), my_rank), std::size_t{0});
    }
};

inline void check_array_name(const std::string& name) {
    if (name.empty())
        throw ValueError("VTKHDF array names must not be empty");
    if (name.find_first_of("/.") != std::string::npos)
        throw ValueError("VTKHDF array names must not contain '/' or '.' (received '" + name + "')");
}

inline bool is_string_field(const Field& field) {
    return field.precision().template is<char>() && field.layout().dimension() == 1;
}

inline std::vector<std::string> extract_strings(const Field& field) {
    const auto strings = field.template export_to<MultiString>();
    const auto slices = strings.slices();
    return {slices.begin(), slices.end()};
}
 
#if GRIDFORMAT_HAVE_HIGH_FIVE

template<typename C>
class DataSetField : public LazyField<const HDF5::File<C>&> {
    using ParentType = LazyField<const HDF5::File<C>&>;
 
 public:
    using ParentType::ParentType;
 
    explicit DataSetField(const HDF5::File<C>& file, std::string path)
    : ParentType{
        file,
        MDLayout{file.get_dimensions(path).value()},
        DynamicPrecision{file.get_precision(path).value()},
        [_p=std::move(path)] (const HDF5::File<C>& file) {
            return file.visit_dataset(_p, [&] <typename F> (F&& field) {
                return field.serialized();
            });
        }
    } {}
};
 
template<typename C>
DataSetField(const HDF5::File<C>&, std::string) -> DataSetField<C>;
 
template<typename C, typename CB>
DataSetField(const HDF5::File<C>&, MDLayout, DynamicPrecision, CB&&) -> DataSetField<C>;

template<typename C>
std::string get_file_type(const HDF5::File<C>& file) {
    if (!file.exists("/VTKHDF"))
        throw IOError("Given file is not a VTK-HDF file");
    if (!file.has_attribute_at("/VTKHDF/Type"))
        throw IOError("VTKHDF-Type attribute missing");
    return file.template read_attribute_to<std::string>("/VTKHDF/Type");
}

template<typename C>
void check_version_compatibility(const HDF5::File<C>& file, const std::array<std::size_t, 2>& supported) {
    if (file.has_attribute_at("/VTKHDF/Version"))
        file.visit_attribute("/VTKHDF/Version", [&] (auto&& field) {
            const auto version = field.template export_to<std::vector<std::size_t>>();
            if (version.size() > 0 && version.at(0) > supported[0])
                throw ValueError(
                    "File version is higher than supported by the reader (" + as_string(supported, ".") + ")"
                );
        });
}

inline FieldPtr as_field_data_tuples(FieldPtr field) {
    const auto layout = field->layout();
    if (layout.dimension() == 0)
        return make_field_ptr(ReshapedField{field, MDLayout{{std::size_t{1}}}});
    if (layout.dimension() > 2)
        return make_field_ptr(ReshapedField{field, MDLayout{{layout.extent(0), layout.number_of_entries(1)}}});
    return field;
}

inline void write_field_data(HDF5::File<>& file, const std::string& name, FieldPtr field) {
    const std::string path = "/VTKHDF/FieldData/" + name;
    if (is_string_field(*field))
        file.write_strings(extract_strings(*field), path);
    else
        file.write(*as_field_data_tuples(field), path);
}

inline void write_field_data_step(HDF5::File<>& file, const std::string& name, FieldPtr field) {
    const std::string path = "/VTKHDF/FieldData/" + name;
    const auto offset = file.get_dimensions(path).value_or(std::vector<std::size_t>{0}).at(0);
    std::array<std::size_t, 2> sizes;
    if (is_string_field(*field)) {
        const auto strings = extract_strings(*field);
        sizes = {1, strings.size()};
        file.write_strings(strings, path);
    } else {
        const auto tuples = as_field_data_tuples(field);
        const auto layout = tuples->layout();
        sizes = {layout.dimension() > 1 ? layout.extent(1) : 1, layout.extent(0)};
        file.write(*tuples, path);
    }
    file.write(std::array{offset}, "/VTKHDF/Steps/FieldDataOffsets/" + name);
    file.write(std::vector{sizes}, "/VTKHDF/Steps/FieldDataSizes/" + name);
}

inline void write_field_data_step_referencing(std::size_t step_index, HDF5::File<>& file, const std::string& name) {
    const std::string offsets_path = "/VTKHDF/Steps/FieldDataOffsets/" + name;
    const std::string sizes_path = "/VTKHDF/Steps/FieldDataSizes/" + name;
    const auto offset = file.read_dataset_to<std::size_t>(
        offsets_path, HDF5::Slice{.offset = {step_index}, .count = {1}}
    );
    const auto sizes = file.read_dataset_to<std::vector<std::size_t>>(
        sizes_path, HDF5::Slice{.offset = {step_index, 0}, .count = {1, 2}}
    );
    file.write(std::array{offset}, offsets_path);
    file.write(std::vector{std::array{sizes.at(0), sizes.at(1)}}, sizes_path);
}

template<typename C>
FieldPtr read_field_data(const HDF5::File<C>& file, const std::string& name, std::optional<std::size_t> step) {
    const std::string path = "VTKHDF/FieldData/" + name;
    const auto precision = file.get_precision(path).value();
    const bool is_string = precision.template is<char>();
    if (!step)
        return is_string ? make_field_ptr(RangeField{file.template read_dataset_to<MultiString>(path)})
                         : make_field_ptr(DataSetField{file, path});
 
    // Without sizes (written by GridFormat <= 0.5), each step is a single row of the dataset
    const auto offset = file.template read_dataset_to<std::size_t>(
        "VTKHDF/Steps/FieldDataOffsets/" + name, HDF5::Slice{.offset = {*step}, .count = {1}}
    );
    const std::string sizes_path = "VTKHDF/Steps/FieldDataSizes/" + name;
    const auto sizes = file.has_dataset_at(sizes_path)  // (components, tuples)
        ? std::optional{file.template read_dataset_to<std::vector<std::size_t>>(
            sizes_path, HDF5::Slice{.offset = {*step, 0}, .count = {1, 2}}
        )}
        : std::nullopt;
    if (is_string) {
        const std::size_t count = sizes ? sizes->at(0)*sizes->at(1) : 1;
        return make_field_ptr(RangeField{file.template read_dataset_to<MultiString>(
            path, HDF5::Slice{.offset = {offset}, .count = {count}}
        )});
    }
 
    const auto dims = file.get_dimensions(path).value();
    std::vector<std::size_t> slice_offset(dims.size(), 0);
    slice_offset.at(0) = offset;
    auto count = dims;
    count.at(0) = sizes ? sizes->at(1) : 1;
    if (sizes && MDLayout{count}.number_of_entries() != sizes->at(0)*sizes->at(1))
        throw SizeError("Field data '" + name + "' does not match its stated number of components");
    MDLayout layout = sizes ? MDLayout{count} : MDLayout{count | std::views::drop(1)};
    return make_field_ptr(DataSetField{
        file,
        std::move(layout),
        precision,
        [p=path, o=std::move(slice_offset), c=std::move(count)] (const HDF5::File<C>& f) {
            return f.visit_dataset(p, [&] <typename F> (F&& field) {
                return FlattenedField{make_field_ptr(std::move(field))}.serialized();
            }, HDF5::Slice{.offset = o, .count = c});
        }
    });
}

template<Concepts::Communicator C, std::invocable<HDF5::File<>&> Write>
void write_on_root(const std::string& filename, const C& comm, int root, Write&& write) {
    Parallel::barrier(comm);
    std::exception_ptr error;
    if (Parallel::rank(comm) == root) {
        try {
            HDF5::File file{filename, HDF5::File<>::append};
            write(file);
        } catch (...) {
            error = std::current_exception();
        }
    }
    // the broadcast also keeps the other ranks from reopening the file before the root is done
    if (Parallel::broadcast(comm, static_cast<int>(error != nullptr), root)) {
        if (error)
            std::rethrow_exception(error);
        throw IOError("Writing to '" + filename + "' failed on rank " + std::to_string(root));
    }
}
 
#endif  // GRIDFORMAT_HAVE_HIGH_FIVE
 
}  // namespace VTKHDF
}  // namespace GridFormat
 
#endif  // GRIDFORMAT_VTK_HDF_COMMON_HPP_