GridFormat 0.4.0
I/O-Library for grid-like data structures
Loading...
Searching...
No Matches
hdf_unstructured_grid_writer.hpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2022-2023 Dennis Gläser <dennis.glaeser@iws.uni-stuttgart.de>
2// SPDX-License-Identifier: MIT
8#ifndef GRIDFORMAT_VTK_HDF_UNSTRUCTURED_GRID_WRITER_HPP_
9#define GRIDFORMAT_VTK_HDF_UNSTRUCTURED_GRID_WRITER_HPP_
10#if GRIDFORMAT_HAVE_HIGH_FIVE
11
12#include <type_traits>
13#include <ostream>
14
15#include <gridformat/common/exceptions.hpp>
16#include <gridformat/common/md_layout.hpp>
17#include <gridformat/common/concepts.hpp>
18#include <gridformat/common/lvalue_reference.hpp>
19
20#include <gridformat/parallel/communication.hpp>
21#include <gridformat/parallel/concepts.hpp>
22
25
28
29namespace GridFormat {
30
31template<bool is_transient,
32 Concepts::UnstructuredGrid G,
33 Concepts::Communicator Communicator = NullCommunicator>
34class VTKHDFUnstructuredGridWriterImpl : public GridDetail::WriterBase<is_transient, G>::type {
35 static constexpr int root_rank = 0;
36 static constexpr std::size_t vtk_space_dim = 3;
37
38 using CT = CoordinateType<G>;
39 template<typename T> using Vector = std::array<T, vtk_space_dim>;
40 template<typename T> using Tensor = std::array<T, vtk_space_dim*vtk_space_dim>;
41
43 using HDF5File = HDF5::File<Communicator>;
44
45 static constexpr WriterOptions writer_opts{
47 .append_null_terminator_to_strings = true
48 };
49
50 struct TimeSeriesOffsets {
51 std::size_t cell_offset;
52 std::size_t connectivity_offset;
53 std::size_t point_offset;
54 };
55
56 public:
57 using Grid = G;
58
59 explicit VTKHDFUnstructuredGridWriterImpl(LValueReferenceOf<const Grid> grid)
60 requires(!is_transient && std::is_same_v<Communicator, NullCommunicator>)
61 : GridWriter<Grid>(grid.get(), ".hdf", writer_opts)
62 {}
63
64 explicit VTKHDFUnstructuredGridWriterImpl(LValueReferenceOf<const Grid> grid, const Communicator& comm)
65 requires(!is_transient && std::is_copy_constructible_v<Communicator>)
66 : GridWriter<Grid>(grid.get(), ".hdf", writer_opts)
67 , _comm{comm}
68 {}
69
70 explicit VTKHDFUnstructuredGridWriterImpl(LValueReferenceOf<const Grid> grid,
71 std::string filename_without_extension,
73 requires(is_transient && std::is_same_v<Communicator, NullCommunicator>)
74 : TimeSeriesGridWriter<Grid>(grid.get(), writer_opts)
75 , _comm{}
76 , _timeseries_filename{std::move(filename_without_extension) + ".hdf"}
77 , _transient_opts{std::move(opts)}
78 {}
79
80 explicit VTKHDFUnstructuredGridWriterImpl(LValueReferenceOf<const Grid> grid,
81 const Communicator& comm,
82 std::string filename_without_extension,
84 requires(is_transient && std::is_copy_constructible_v<Communicator>)
85 : TimeSeriesGridWriter<Grid>(grid.get(), writer_opts)
86 , _comm{comm}
87 , _timeseries_filename{std::move(filename_without_extension) + ".hdf"}
88 , _transient_opts{std::move(opts)}
89 {}
90
91 const Communicator& communicator() const {
92 return _comm;
93 }
94
95 private:
96 void _write(std::ostream&) const {
97 throw InvalidState("VTKHDFUnstructuredGridWriter does not support export into stream");
98 }
99
100 void _write(const std::string& filename_with_ext) const {
101 if constexpr (is_transient)
102 throw InvalidState("This overload only works for non-transient output");
103 HDF5File file{filename_with_ext, _comm, HDF5File::overwrite};
104 _write_to(file);
105 }
106
107 std::string _write(double t) {
108 if constexpr (!is_transient)
109 throw InvalidState("This overload only works for transient output");
110
111 if (this->_step_count == 0)
112 HDF5File::clear(_timeseries_filename, _comm);
113
114 HDF5File file{_timeseries_filename, _comm, HDF5File::append};
115 const auto offsets = _write_to(file);
116
117 file.write_attribute(this->_step_count+1, "/VTKHDF/Steps/NSteps");
118 file.write(std::array{t}, "VTKHDF/Steps/Values");
119 file.write(std::vector{offsets.point_offset}, "VTKHDF/Steps/PointOffsets");
120 file.write(std::vector{std::array{offsets.cell_offset}}, "/VTKHDF/Steps/CellOffsets");
121 file.write(std::vector{std::array{offsets.connectivity_offset}}, "VTKHDF/Steps/ConnectivityIdOffsets");
122
123 file.write(std::vector{Parallel::size(_comm)}, "/VTKHDF/Steps/NumberOfParts");
124 if (this->_step_count > 0 && _transient_opts.static_grid) {
125 file.write(
126 std::vector{_get_last_step_data(file, "PartOffsets")},
127 "/VTKHDF/Steps/PartOffsets"
128 );
129 } else {
130 const std::size_t offset = this->_step_count == 0
131 ? 0
132 : _get_last_step_data(file, "PartOffsets") + Parallel::size(_comm);
133 file.write(std::vector{offset}, "/VTKHDF/Steps/PartOffsets");
134 }
135
136 return _timeseries_filename;
137 }
138
139 TimeSeriesOffsets _write_to(HDF5File& file) const {
140 file.write_attribute(std::array<std::size_t, 2>{(is_transient ? 2 : 1), 0}, "/VTKHDF/Version");
141 file.write_attribute("UnstructuredGrid", "/VTKHDF/Type");
142
143 TimeSeriesOffsets offsets;
144 const auto context = IOContext::from(this->grid(), _comm, root_rank);
145 _write_num_cells_and_points(file, context);
146 offsets.point_offset = _write_coordinates(file, context);
147 offsets.connectivity_offset = _write_connectivity(file, context);
148 offsets.cell_offset = _write_types(file, context);
149 _write_offsets(file, context);
150 _write_meta_data(file);
151 _write_point_fields(file, context);
152 _write_cell_fields(file, context);
153
154 return offsets;
155 }
156
157 void _write_num_cells_and_points(HDF5File& file, const IOContext& context) const {
158 _write_values(file, "/VTKHDF/NumberOfPoints", std::vector{number_of_points(this->grid())}, context);
159 _write_values(file, "/VTKHDF/NumberOfCells", std::vector{number_of_cells(this->grid())}, context);
160 }
161
162 std::size_t _write_coordinates(HDF5File& file, const IOContext& context) const {
163 if constexpr (is_transient) {
164 if (this->_step_count > 0 && _transient_opts.static_grid)
165 return _get_last_step_data(file, "PointOffsets");
166 }
167 const auto coords_field = VTK::make_coordinates_field<CT>(this->grid(), false);
168 const auto offset = _get_current_offset(file, "/VTKHDF/Points");
169 _write_point_field(file, "/VTKHDF/Points", *coords_field, context);
170 return offset;
171 }
172
173 std::size_t _write_connectivity(HDF5File& file, const IOContext& context) const {
174 if constexpr (is_transient) {
175 if (this->_step_count > 0 && _transient_opts.static_grid)
176 return _get_last_step_data(file, "ConnectivityIdOffsets");
177 }
178 const auto point_id_map = make_point_id_map(this->grid());
179 const auto connectivity_field = VTK::make_connectivity_field(this->grid(), point_id_map);
180 const auto num_entries = connectivity_field->layout().number_of_entries();
181 const auto my_num_ids = std::vector{static_cast<long>(num_entries)};
182 std::vector<long> connectivity(num_entries);
183 connectivity_field->export_to(connectivity);
184 const auto offset = _get_current_offset(file, "/VTKHDF/Connectivity");
185 _write_values(file, "/VTKHDF/Connectivity", connectivity, context);
186 _write_values(file, "/VTKHDF/NumberOfConnectivityIds", my_num_ids, context);
187 return offset;
188 }
189
190 std::size_t _write_types(HDF5File& file, const IOContext& context) const {
191 if constexpr (is_transient) {
192 if (this->_step_count > 0 && _transient_opts.static_grid)
193 return _get_last_step_data(file, "CellOffsets");
194 }
195 const auto types_field = VTK::make_cell_types_field(this->grid());
196 std::vector<std::uint8_t> types(types_field->layout().number_of_entries());
197 types_field->export_to(types);
198 const auto offset = _get_current_offset(file, "VTKHDF/Types");
199 _write_values(file, "VTKHDF/Types", types, context);
200 return offset;
201 }
202
203 std::size_t _write_offsets(HDF5File& file, const IOContext& context) const {
204 if constexpr (is_transient) {
205 if (this->_step_count > 0 && _transient_opts.static_grid)
206 return _get_last_step_data(file, "CellOffsets");
207 }
208 const auto offsets_field = VTK::make_offsets_field(this->grid());
209 const auto num_offset_entries = offsets_field->layout().number_of_entries() + 1;
210 std::vector<long> offsets(num_offset_entries);
211 offsets_field->export_to(std::ranges::subrange(std::next(offsets.begin()), offsets.end()));
212 offsets[0] = long{0};
213 const auto offset = _get_current_offset(file, "/VTKHDF/Offsets");
214 _write_values(file, "/VTKHDF/Offsets", offsets, context);
215 return offset;
216 }
217
218 void _write_meta_data(HDF5File& file) const {
219 std::ranges::for_each(this->_meta_data_field_names(), [&] (const std::string& name) {
220 if constexpr (is_transient) {
221 if (this->_step_count > 0 && _transient_opts.static_meta_data) {
222 file.write(std::array{0}, "/VTKHDF/Steps/FieldDataOffsets/" + name);
223 return;
224 } else {
225 file.write(std::array{this->_step_count}, "/VTKHDF/Steps/FieldDataOffsets/" + name);
226 }
227 // For transient data, prepend a dimension indicating the step count
228 TransformedField sub{this->_get_meta_data_field_ptr(name), FieldTransformation::as_sub_field};
229 file.write(sub, "/VTKHDF/FieldData/" + name);
230 } else {
231 file.write(*this->_get_meta_data_field_ptr(name), "/VTKHDF/FieldData/" + name);
232 }
233 });
234 }
235
236 void _write_point_fields(HDF5File& file, const IOContext& context) const {
237 std::ranges::for_each(this->_point_field_names(), [&] (const std::string& name) {
238 if constexpr (is_transient)
239 _write_step_offset(
240 file,
241 _get_current_offset(file, "/VTKHDF/PointData/" + name),
242 "/VTKHDF/Steps/PointDataOffsets/" + name
243 );
244 auto reshaped = _reshape(VTK::make_vtk_field(this->_get_point_field_ptr(name)));
245 _write_point_field(file, "/VTKHDF/PointData/" + name, *reshaped, context);
246 });
247 }
248
249 void _write_cell_fields(HDF5File& file, const IOContext& context) const {
250 std::ranges::for_each(this->_cell_field_names(), [&] (const std::string& name) {
251 if constexpr (is_transient)
252 _write_step_offset(
253 file,
254 _get_current_offset(file, "VTKHDF/CellData/" + name),
255 "/VTKHDF/Steps/CellDataOffsets/" + name
256 );
257 auto reshaped = _reshape(VTK::make_vtk_field(this->_get_cell_field_ptr(name)));
258 _write_cell_field(file, "/VTKHDF/CellData/" + name, *reshaped, context);
259 });
260 }
261
262 FieldPtr _reshape(FieldPtr field_ptr) const {
263 // VTK requires tensors to be written as flat fields
264 const auto layout = field_ptr->layout();
265 if (layout.dimension() > 2)
266 return make_field_ptr(ReshapedField{
267 field_ptr,
268 MDLayout{{layout.extent(0), layout.number_of_entries(1)}}
269 });
270 return field_ptr;
271 }
272
273 template<Concepts::Scalar T>
274 void _write_values(HDF5File& file,
275 const std::string& path,
276 const std::vector<T>& values,
277 const IOContext& context) const {
278 if (context.is_parallel) {
279 const auto num_values = values.size();
280 const auto total_num_values = Parallel::sum(_comm, num_values, root_rank);
281 const auto my_total_num_values = Parallel::broadcast(_comm, total_num_values, root_rank);
282 const auto my_offset = _accumulate_rank_offset(num_values);
283 HDF5::Slice slice{
284 .offset = std::vector{my_offset},
285 .count = std::vector{num_values},
286 .total_size = std::vector{my_total_num_values}
287 };
288 file.write(values, path, slice);
289 } else {
290 file.write(values, path);
291 }
292 }
293
294 void _write_point_field(HDF5File& file,
295 const std::string& path,
296 const Field& field,
297 const IOContext& context) const {
298 _write_field(file, path, field, context.is_parallel, context.my_point_offset, context.num_points_total);
299 }
300
301 void _write_cell_field(HDF5File& file,
302 const std::string& path,
303 const Field& field,
304 const IOContext& context) const {
305 _write_field(file, path, field, context.is_parallel, context.my_cell_offset, context.num_cells_total);
306 }
307
308 void _write_field(HDF5File& file,
309 const std::string& path,
310 const Field& field,
311 bool is_parallel,
312 std::size_t main_offset,
313 std::size_t main_size) const {
314 if (is_parallel) {
315 const auto layout = field.layout();
316 std::vector<std::size_t> count(layout.dimension());
317 layout.export_to(count);
318
319 std::vector<std::size_t> size = count;
320 std::vector<std::size_t> offset(layout.dimension(), 0);
321 offset.at(0) = main_offset;
322 size.at(0) = main_size;
323
324 file.write(field, path, HDF5::Slice{
325 .offset = std::move(offset),
326 .count = std::move(count),
327 .total_size = std::move(size)
328 });
329 } else {
330 file.write(field, path);
331 }
332 }
333
334 std::size_t _accumulate_rank_offset(std::size_t my_size) const {
335 auto all_offsets = Parallel::gather(_comm, my_size, root_rank);
336 if (Parallel::rank(_comm) == root_rank) {
337 std::partial_sum(all_offsets.begin(), std::prev(all_offsets.end()), all_offsets.begin());
338 std::copy(std::next(all_offsets.rbegin()), all_offsets.rend(), all_offsets.rbegin());
339 all_offsets[0] = 0;
340 }
341 const auto my_offset = Parallel::scatter(_comm, all_offsets, root_rank);
342 if (my_offset.size() != 1)
343 throw ValueError("Unexpected scatter result");
344 return my_offset[0];
345 }
346
347 void _write_step_offset(HDF5File& file,
348 const std::size_t offset,
349 const std::string& path) const {
350 file.write(std::array{offset}, path);
351 }
352
353 std::size_t _get_current_offset(HDF5File& file, const std::string& path) const {
354 const auto cur_dimensions = file.get_dimensions(path);
355 return cur_dimensions ? (*cur_dimensions)[0] : 0;
356 }
357
358 std::size_t _get_last_step_data(const HDF5File& file, const std::string& sub_path) const {
359 if (this->_step_count == 0)
360 throw ValueError("Last step data can only be read after at least one write");
361 const std::string path ="VTKHDF/Steps/" + sub_path;
362
363 auto step_dimensions = file.get_dimensions(path).value();
364 std::ranges::fill(step_dimensions, std::size_t{1});
365
366 auto access_offset = step_dimensions;
367 std::ranges::fill(access_offset, std::size_t{0});
368 access_offset.at(0) = this->_step_count - 1;
369
370 return file.template read_dataset_to<std::size_t>(path, HDF5::Slice{
371 .offset = std::move(access_offset),
372 .count = std::move(step_dimensions)
373 });
374 }
375
376 Communicator _comm;
377 std::string _timeseries_filename = "";
378 VTK::HDFTransientOptions _transient_opts;
379};
380
385template<Concepts::UnstructuredGrid G, Concepts::Communicator C = NullCommunicator>
388 public:
389 using ParentType::ParentType;
390};
391
396template<Concepts::UnstructuredGrid G, Concepts::Communicator C = NullCommunicator>
399 public:
400 using ParentType::ParentType;
401};
402
403template<Concepts::UnstructuredGrid G>
405template<Concepts::UnstructuredGrid G, Concepts::Communicator C>
407
408template<Concepts::UnstructuredGrid G>
409VTKHDFUnstructuredTimeSeriesWriter(const G&, std::string, VTK::HDFTransientOptions = {}) -> VTKHDFUnstructuredTimeSeriesWriter<G, NullCommunicator>;
410template<Concepts::UnstructuredGrid G, Concepts::Communicator C>
411VTKHDFUnstructuredTimeSeriesWriter(const G&, const C&, std::string, VTK::HDFTransientOptions = {}) -> VTKHDFUnstructuredTimeSeriesWriter<G, C>;
412
413} // namespace GridFormat
414
415#endif // GRIDFORMAT_HAVE_HIGH_FIVE
416#endif // GRIDFORMAT_VTK_HDF_UNSTRUCTURED_GRID_WRITER_HPP_
Abstract interface for fields of values that is used by writers/readers to store fields.
Definition: field.hpp:37
MDLayout layout() const
Return the layout of this field.
Definition: field.hpp:51
Abstract base class for grid file writers.
Definition: writer.hpp:303
Abstract base class for time series file writers.
Definition: writer.hpp:344
Definition: hdf_unstructured_grid_writer.hpp:34
Writer for the VTK HDF file format for unstructured grids.
Definition: hdf_unstructured_grid_writer.hpp:386
Writer for the transient VTK HDF file format for unstructured grids.
Definition: hdf_unstructured_grid_writer.hpp:397
Grid concepts.
Base classes for grid data writers.
std::shared_ptr< const Field > FieldPtr
Pointer type used by writers/readers for fields.
Definition: field.hpp:186
FieldPtr make_field_ptr(F &&f)
Factory function for field pointers.
Definition: field.hpp:192
Common functionality for writing VTK HDF files.
Helper class to store processor offsets in parallel writes.
Definition: hdf_common.hpp:37
Options for transient vtk-hdf file formats.
Definition: hdf_common.hpp:26
bool static_grid
Set to true the grid is the same for all time steps (will only be written once)
Definition: hdf_common.hpp:27
bool static_meta_data
Set to true if the metadata is same for all time steps (will only be written once)
Definition: hdf_common.hpp:28
Options that writer implementations can pass to the base class.
Definition: writer.hpp:59
bool use_structured_grid_ordering
Use row-major structured grid ordering.
Definition: writer.hpp:60
Common functionality for VTK writers.