kspaceFirstOrder3D-CUDA  1.1
The CUDA/C++ implementation of the k-wave toolbox for the time-domain simulation of acoustic wave fields in 3D
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CuboidOutputHDF5Stream.cpp
Go to the documentation of this file.
1 /**
2  * @file CuboidOutputHDF5Stream.cpp
3  *
4  * @author Jiri Jaros \n
5  * Faculty of Information Technology \n
6  * Brno University of Technology \n
7  * jarosjir@fit.vutbr.cz
8  *
9  * @brief The implementation file of classes responsible for storing output quantities based
10  * on the cuboid sensor mask into the output HDF5 file.
11  *
12  * @version kspaceFirstOrder3D 3.4
13  *
14  * @date 13 February 2015, 12:51 (created)
15  * 29 July 2016, 16:56 (revised)
16  *
17  * @section License
18  * This file is part of the C++ extension of the k-Wave Toolbox
19  * (http://www.k-wave.org).\n Copyright (C) 2016 Jiri Jaros and Bradley Treeby.
20  *
21  * This file is part of the k-Wave. k-Wave is free software: you can redistribute it and/or modify
22  * it under the terms of the GNU Lesser General Public License as published by the Free Software
23  * Foundation, either version 3 of the License, or (at your option) any later version.
24  *
25  * k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even
26  * the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
27  * General Public License for more details.
28  *
29  * You should have received a copy of the GNU Lesser General Public License along with k-Wave.
30  * If not, see http://www.gnu.org/licenses/.
31  */
32 
35 
36 #include <Parameters/Parameters.h>
37 #include <Logger/Logger.h>
38 
39 using std::string;
40 
41 //------------------------------------------------------------------------------------------------//
42 //------------------------------------------ Constants -------------------------------------------//
43 //------------------------------------------------------------------------------------------------//
44 
45 //------------------------------------------------------------------------------------------------//
46 //--------------------------------------- Public methods -----------------------------------------//
47 //------------------------------------------------------------------------------------------------//
48 
49 /**
50  * Constructor - links the HDF5 dataset, source (sampled matrix), sensor mask and the reduce
51  * operator together. The constructor DOES NOT allocate memory because the size of the sensor mask
52  * is not known at the time the instance of the class is being created.
53  *
54  * @param [in] file - HDF5 file to write the output to
55  * @param [in] groupName - The name of the HDF5 group with datasets for particular cuboids
56  * @param [in] sourceMatrix - Source matrix to be sampled
57  * @param [in] sensorMask - Sensor mask with the cuboid coordinates
58  * @param [in] reduceOp - Reduce operator
59 
60  */
62  TMatrixName& groupName,
63  const TRealMatrix& sourceMatrix,
64  const TIndexMatrix& sensorMask,
65  const TReduceOperator reduceOp)
66  : TBaseOutputHDF5Stream(file, groupName, sourceMatrix, reduceOp),
67  sensorMask(sensorMask),
68  group(H5I_BADID),
69  sampledTimeStep(0),
70  eventSamplingFinished()
71 {
72  // Create event for sampling
73  checkCudaErrors(cudaEventCreate(&eventSamplingFinished));
74 }// end of TCubodidOutputHDF5Stream
75 //--------------------------------------------------------------------------------------------------
76 
77 /**
78  * Destructor.
79  * if the file is still opened, it applies the post processing and flush the data.
80  * Then, the object memory is freed and the object destroyed.
81  */
83 {
84  // Destroy sampling event
85  checkCudaErrors(cudaEventDestroy(eventSamplingFinished));
86  // Close the stream
87  Close();
88  // free memory
89  FreeMemory();
90 }// end ~TCubodidOutputHDF5Stream
91 //--------------------------------------------------------------------------------------------------
92 
93 /**
94  * Create a HDF5 stream and allocate data for it. It also creates a HDF5 group with particular
95  * datasets (one per cuboid).
96  */
98 {
99  // Create the HDF5 group and open it
101 
102  // Create all datasets (sizes, chunks, and attributes)
103  size_t nCuboids = sensorMask.GetDimensionSizes().ny;
104  cuboidsInfo.reserve(nCuboids);
105  size_t actualPositionInBuffer = 0;
106 
107  for (size_t cuboidIdx = 0; cuboidIdx < nCuboids; cuboidIdx++)
108  {
109  TCuboidInfo cuboidInfo;
110 
111  cuboidInfo.cuboidIdx = CreateCuboidDataset(cuboidIdx);
112  cuboidInfo.startingPossitionInBuffer = actualPositionInBuffer;
113  cuboidsInfo.push_back(cuboidInfo);
114 
115  actualPositionInBuffer += (sensorMask.GetBottomRightCorner(cuboidIdx) -
116  sensorMask.GetTopLeftCorner(cuboidIdx)
117  ).GetElementCount();
118  }
119 
120  //we're at the beginning
121  sampledTimeStep = 0;
122 
123  // Create the memory buffer if necessary and set starting address
125 
126  // Allocate memory
127  AllocateMemory();
128 }// end of Create
129 //--------------------------------------------------------------------------------------------------
130 
131 
132 
133 /**
134  * Reopen the output stream after restart and reload data.
135  */
137 {
138  // Get parameters
139  const TParameters& params = TParameters::GetInstance();
140 
141  sampledTimeStep = 0;
142  if (reduceOp == NONE) // set correct sampled timestep for raw data series
143  {
144  sampledTimeStep = (params.Get_t_index() < params.GetStartTimeIndex()) ?
145  0 : (params.Get_t_index() - params.GetStartTimeIndex());
146  }
147 
148  // Create the memory buffer if necessary and set starting address
150 
151  // Allocate memory if needed
152  AllocateMemory();
153 
154  // Open all datasets (sizes, chunks, and attributes)
155  size_t nCuboids = sensorMask.GetDimensionSizes().ny;
156  cuboidsInfo.reserve(nCuboids);
157  size_t actualPositionInBuffer = 0;
158 
159  // Open the HDF5 group
161 
162  for (size_t CuboidIdx = 0; CuboidIdx < nCuboids; CuboidIdx++)
163  {
164  TCuboidInfo cuboidInfo;
165  // Indexed from 1
166  const string datasetName = std::to_string(CuboidIdx + 1);
167 
168  // open the dataset
169  cuboidInfo.cuboidIdx = file.OpenDataset(group,datasetName.c_str());
170  cuboidInfo.startingPossitionInBuffer = actualPositionInBuffer;
171  cuboidsInfo.push_back(cuboidInfo);
172 
173  // read only if there is anything to read
174  if (params.Get_t_index() > params.GetStartTimeIndex())
175  {
176  if (reduceOp != NONE)
177  { // Reload data
178  TDimensionSizes cuboidSize((sensorMask.GetBottomRightCorner(CuboidIdx) - sensorMask.GetTopLeftCorner(CuboidIdx)).nx,
179  (sensorMask.GetBottomRightCorner(CuboidIdx) - sensorMask.GetTopLeftCorner(CuboidIdx)).ny,
180  (sensorMask.GetBottomRightCorner(CuboidIdx) - sensorMask.GetTopLeftCorner(CuboidIdx)).nz);
181 
183  datasetName.c_str(),
184  cuboidSize,
185  hostBuffer + actualPositionInBuffer);
186  }
187  }
188  // move the pointer for the next cuboid beginning (this inits the locations)
189  actualPositionInBuffer += (sensorMask.GetBottomRightCorner(CuboidIdx) -
190  sensorMask.GetTopLeftCorner(CuboidIdx)).GetElementCount();
191  }
192 
193  // copy data over to the GPU only if there is anything to read
194  if (params.Get_t_index() > params.GetStartTimeIndex())
195  {
197  }
198 }// end of Reopen
199 //--------------------------------------------------------------------------------------------------
200 
201 
202 /**
203  * Sample grid points, line them up in the buffer, if necessary a reduce operator is applied.
204  *
205  * @warning data is not flushed, there is no sync.
206  */
208 {
209  size_t cuboidInBufferStart = 0;
210  // dimension sizes of the matrix being sampled
211  const dim3 dimSizes (sourceMatrix.GetDimensionSizes().nx,
214 
215  // Run over all cuboids - this is not a good solution as we need to run a distinct kernel for a cuboid
216  for (size_t cuboidIdx = 0; cuboidIdx < cuboidsInfo.size(); cuboidIdx++)
217  {
218  // copy down dim sizes
219  const dim3 topLeftCorner(sensorMask.GetTopLeftCorner(cuboidIdx).nx,
220  sensorMask.GetTopLeftCorner(cuboidIdx).ny,
221  sensorMask.GetTopLeftCorner(cuboidIdx).nz);
222  const dim3 bottomRightCorner(sensorMask.GetBottomRightCorner(cuboidIdx).nx,
224  sensorMask.GetBottomRightCorner(cuboidIdx).nz);
225 
226  //get number of samples within the cuboid
227  const size_t nSamples = (sensorMask.GetBottomRightCorner(cuboidIdx) -
228  sensorMask.GetTopLeftCorner(cuboidIdx)
229  ).GetElementCount();
230 
231  switch (reduceOp)
232  {
233  case NONE :
234  {
235  // Kernel to sample raw quantities inside one cuboid
236  OutputStreamsCUDAKernels::SampleCuboid<NONE>
237  (deviceBuffer + cuboidInBufferStart,
239  topLeftCorner,
240  bottomRightCorner,
241  dimSizes,
242  nSamples);
243  break;
244  }
245  case RMS :
246  {
247  OutputStreamsCUDAKernels::SampleCuboid<RMS>
248  (deviceBuffer + cuboidInBufferStart,
250  topLeftCorner,
251  bottomRightCorner,
252  dimSizes,
253  nSamples);
254  break;
255  }
256  case MAX :
257  {
258  OutputStreamsCUDAKernels::SampleCuboid<MAX>
259  (deviceBuffer + cuboidInBufferStart,
261  topLeftCorner,
262  bottomRightCorner,
263  dimSizes,
264  nSamples);
265  break;
266  }
267  case MIN :
268  {
269  OutputStreamsCUDAKernels::SampleCuboid<MIN>
270  (deviceBuffer + cuboidInBufferStart,
272  topLeftCorner,
273  bottomRightCorner,
274  dimSizes,
275  nSamples);
276  break;
277  }
278  }
279 
280  cuboidInBufferStart += nSamples;
281  }
282 
283  if (reduceOp == NONE)
284  {
285  // Record an event when the data has been copied over.
286  checkCudaErrors(cudaEventRecord(eventSamplingFinished));
287  }
288 }// end of Sample
289 //-------------------------------------------------------------------------------------------------
290 
291 /**
292  * Flush data for the timestep. Only applicable on RAW data series.
293  */
295 {
296  if (reduceOp == NONE)
297  {
298  // make sure the data has been copied from the GPU
299  cudaEventSynchronize(eventSamplingFinished);
300 
301  // only raw time series are flushed down to the disk every time step
303  }
304 }// end of FlushRaw
305 //--------------------------------------------------------------------------------------------------
306 
307 
308 /*
309  * Apply post-processing on the buffer and flush it to the file.
310  */
312 {
313  // run inherited method
315 
316  // When no reduce operator is applied, the data is flushed after every time step
317  // which means it has been done before
318  if (reduceOp != NONE)
319  {
320  // Copy data from GPU matrix
322 
324  }
325 }// end of PostProcessing
326 //--------------------------------------------------------------------------------------------------
327 
328 
329 /**
330  * Checkpoint the stream and close.
331  */
333 {
334  // raw data has already been flushed, others has to be flushed here
335  if (reduceOp != NONE)
336  {
337  // copy data from the device
339  // flush to disk
341  }
342 }// end of Checkpoint
343 //-------------------------------------------------------------------------------------------------
344 
345 
346 /**
347  * Close stream (apply post-processing if necessary, flush data, close datasets and the group).
348  */
350 {
351  // the group is still open
352  if (group != H5I_BADID)
353  {
354  // Close all datasets and the group
355  for (size_t cuboidIdx = 0; cuboidIdx < cuboidsInfo.size(); cuboidIdx++)
356  {
357  file.CloseDataset(cuboidsInfo[cuboidIdx].cuboidIdx);
358  }
359  cuboidsInfo.clear();
360 
362  group = H5I_BADID;
363  }// if opened
364 }// end of Close
365 //--------------------------------------------------------------------------------------------------
366 
367 
368 //------------------------------------------------------------------------------------------------//
369 //-------------------------------------- Protected methods ---------------------------------------//
370 //------------------------------------------------------------------------------------------------//
371 
372 /**
373  * Create a new dataset for a given cuboid specified by index (order).
374  *
375  * @param [in] cuboidIdx - Index of the cuboid in the sensor mask
376  * @return HDF5 handle to the dataset.
377  */
378 hid_t TCuboidOutputHDF5Stream::CreateCuboidDataset(const size_t cuboidIdx)
379 {
380  const TParameters& params = TParameters::GetInstance();
381 
382  // if time series then Number of steps else 1
383  const size_t nSampledTimeSteps = (reduceOp == NONE)
384  ? params.Get_nt() - params.GetStartTimeIndex() : 0; // will be a 3D dataset
385 
386  // Set cuboid dimensions (subtract two corners (add 1) and use the appropriate component)
387  TDimensionSizes cuboidSize((sensorMask.GetBottomRightCorner(cuboidIdx) - sensorMask.GetTopLeftCorner(cuboidIdx)).nx,
388  (sensorMask.GetBottomRightCorner(cuboidIdx) - sensorMask.GetTopLeftCorner(cuboidIdx)).ny,
389  (sensorMask.GetBottomRightCorner(cuboidIdx) - sensorMask.GetTopLeftCorner(cuboidIdx)).nz,
390  nSampledTimeSteps);
391 
392  // Set chunk size
393  // If the size of the cuboid is bigger than 32 MB per timestep, set the chunk to approx 4MB
394  size_t nSlabs = 1; //at least one slab
395  TDimensionSizes cuboidChunkSize(cuboidSize.nx, cuboidSize.ny, cuboidSize.nz, (reduceOp == NONE) ? 1 : 0);
396 
397  if (cuboidChunkSize.GetElementCount() > (CHUNK_SIZE_4MB * 8))
398  {
399  while (nSlabs * cuboidSize.nx * cuboidSize.ny < CHUNK_SIZE_4MB) nSlabs++;
400  cuboidChunkSize.nz = nSlabs;
401  }
402 
403  // Indexed from 1
404  const string datasetName = std::to_string(cuboidIdx + 1);
405 
407  datasetName.c_str(),
408  cuboidSize,
409  cuboidChunkSize,
410  params.GetCompressionLevel());
411 
412  // Write dataset parameters
413  file.WriteMatrixDomainType(group, datasetName.c_str(), THDF5_File::REAL);
414  file.WriteMatrixDataType (group, datasetName.c_str(), THDF5_File::FLOAT);
415 
416  return dataset;
417 }//end of CreateCuboidDatasets
418 //--------------------------------------------------------------------------------------------------
419 
420 /**
421  * Flush the buffer to the file (to multiple datasets if necessary).
422  */
424 {
425  TDimensionSizes position (0,0,0,0);
426  TDimensionSizes blockSize(0,0,0,0);
427 
428  if (reduceOp == NONE) position.nt = sampledTimeStep;
429 
430  for (size_t cuboidIdx = 0; cuboidIdx < cuboidsInfo.size(); cuboidIdx++)
431  {
432  blockSize = sensorMask.GetBottomRightCorner(cuboidIdx) - sensorMask.GetTopLeftCorner(cuboidIdx);
433  blockSize.nt = 1;
434 
435  file.WriteHyperSlab(cuboidsInfo[cuboidIdx].cuboidIdx,
436  position,
437  blockSize,
438  hostBuffer + cuboidsInfo[cuboidIdx].startingPossitionInBuffer);
439  }
440 
441  sampledTimeStep++;
442 }// end of FlushBufferToFile
443 //--------------------------------------------------------------------------------------------------
444 
445 //------------------------------------------------------------------------------------------------//
446 //--------------------------------------- Private methods ----------------------------------------//
447 //------------------------------------------------------------------------------------------------//
size_t nx
number of elements in the x direction
size_t GetTotalNumberOfElementsInAllCuboids() const
Get the total number of elements to be sampled within all cuboids.
THDF5_File & file
HDF5 file handle.
hid_t CreateFloatDataset(const hid_t parentGroup, TMatrixName &datasetName, const TDimensionSizes &dimensionSizes, const TDimensionSizes &chunkSizes, const size_t compressionLevel)
Create a float HDF5 dataset at a specified place in the file tree (3D/4D).
Definition: HDF5_File.cpp:296
hid_t group
Handle to a HDF5 dataset.
size_t GetStartTimeIndex() const
Get start time index for sensor recording.
Definition: Parameters.h:225
void WriteMatrixDomainType(const hid_t parentGroup, TMatrixName &datasetName, const TMatrixDomainType &matrixDomainType)
Write matrix domain type into the dataset under the root group.
Definition: HDF5_File.cpp:1048
size_t startingPossitionInBuffer
Having a single buffer for all cuboids, where this one starts.
size_t nt
Number of time steps (for time series datasets).
TDimensionSizes GetBottomRightCorner(const size_t &index) const
Get the bottom right corner of the index-th cuboid.
virtual void CopyDataFromDevice()
Copy data deviceBuffer -> hostBuffer.
calculate root mean square
std::string rootObjectName
Dataset name.
static TParameters & GetInstance()
Get instance of the singleton class.
Definition: Parameters.cpp:70
virtual struct TDimensionSizes GetDimensionSizes() const
Get dimension sizes of the matrix.
hid_t CreateGroup(const hid_t parentGroup, TMatrixName &groupName)
Create a HDF5 group at a specified place in the file tree.
Definition: HDF5_File.cpp:208
const TReduceOperator reduceOp
Reduce operator.
The header file containing the parameters of the simulation.
The header file of cuda kernels used for data sampling (output streams).
static const size_t CHUNK_SIZE_4MB
chunk size of 4MB in number of float elements
const std::string TMatrixName
Datatype for matrix names.
Definition: MatrixNames.h:45
virtual hid_t CreateCuboidDataset(const size_t cuboidIdx)
Create a new dataset for a given cuboid specified by index (order).
#define checkCudaErrors(val)
Macro checking cuda errors and printing the file name and line. Inspired by CUDA common checking rout...
Definition: Logger.h:209
virtual void PostProcess()
Apply post-processing on the buffer and flush it to the file.
hid_t dataset
HDF5 dataset handle.
virtual void AllocateMemory()
A generic function to allocate memory - not used in the base class.
virtual TDimensionSizes GetDimensionSizes() const
Get dimension sizes of the matrix.
void WriteMatrixDataType(const hid_t parentGroup, TMatrixName &datasetName, const TMatrixDataType &matrixDataType)
Write matrix data type into the dataset under a specified group.
Definition: HDF5_File.cpp:1029
const TRealMatrix & sourceMatrix
Source matrix to be sampled.
size_t Get_t_index() const
Get simulation time step.
Definition: Parameters.h:102
The header file containing a class responsible for printing out info and error messages (stdout...
Class storing all parameters of the simulation.
Definition: Parameters.h:49
TDimensionSizes GetTopLeftCorner(const size_t &index) const
Get the top left corner of the index-th cuboid.
virtual void PostProcess()
Apply post-processing on the buffer and flush it to the file.
virtual void FlushRaw()
Flush data to disk (from raw streams only).
hid_t GetRootGroup() const
Get handle to the root group.
Definition: HDF5_File.h:573
std::vector< TCuboidInfo > cuboidsInfo
vector keeping handles and positions of all cuboids
virtual void Create()
Create a HDF5 stream and allocate data for it.
TReduceOperator
How to aggregate data.
float * deviceBuffer
Temporary buffer on the GPU side - only for aggregated quantities.
store actual data (time series)
size_t GetElementCount() const
Get element count, in 3D only spatial domain, in 4D with time.
size_t ny
number of elements in the y direction
cudaEvent_t eventSamplingFinished
Has the sampling finished?
Abstract base class for output data streams (sampled data).
virtual void Checkpoint()
Checkpoint the stream and close.
virtual void FreeMemory()
A generic function to free memory - not used in the base class.
size_t sampledTimeStep
Timestep to store (N/A for aggregated).
The class for real matrices.
Definition: RealMatrix.h:45
float * hostBuffer
Temporary buffer for store on the GPU side.
This structure information about one cuboid. Namely, its HDF5_ID, starting position in a lineup buffe...
virtual void CopyDataToDevice()
Copy data hostBuffer-> deviceBuffer.
const TIndexMatrix & sensorMask
Sensor mask to sample data.
The class for 64b unsigned integers (indices). It is used for sensor_mask_index or sensor_corners_mas...
Definition: IndexMatrix.h:46
TDimensionSizes position
Position in the dataset.
virtual void Sample()
Sample data into buffer, apply reduce or copy to the CPU side.
virtual void Close()
Close stream (apply post-processing if necessary, flush data and close).
size_t GetCompressionLevel() const
Get compression level.
Definition: Parameters.h:218
void WriteHyperSlab(const hid_t dataset, const TDimensionSizes &position, const TDimensionSizes &size, const float *data)
Write a hyper-slab into the dataset - float dataset.
Definition: HDF5_File.cpp:473
virtual float * GetDeviceData()
Get raw GPU data out of the class (for direct GPU kernel access).
The header file of classes responsible for storing output quantities based on the cuboid sensor mask ...
size_t bufferSize
Buffer size.
virtual void Reopen()
Reopen the output stream after restart and reload data.
hid_t OpenGroup(const hid_t parentGroup, TMatrixName &groupName)
Open a HDF5 group at a specified place in the file tree.
Definition: HDF5_File.cpp:232
size_t nz
number of elements in the z direction
void CloseGroup(const hid_t group)
Close group.
Definition: HDF5_File.cpp:254
void ReadCompleteDataset(const hid_t parentGroup, TMatrixName &datasetName, const TDimensionSizes &dimensionSizes, float *data)
Read data from the dataset under a specified group, float dataset.
Definition: HDF5_File.cpp:887
size_t Get_nt() const
Get Nt value.
Definition: Parameters.h:100
virtual void FlushBufferToFile()
Flush the buffer to the file.
void CloseDataset(const hid_t dataset)
Close the HDF5 dataset.
Definition: HDF5_File.cpp:457
TCuboidOutputHDF5Stream(THDF5_File &file, TMatrixName &groupName, const TRealMatrix &sourceMatrix, const TIndexMatrix &sensorMask, const TReduceOperator reduceOp)
Constructor.
hid_t cuboidIdx
Idx of the dataset storing the given cuboid.
Class wrapping the HDF5 routines.
Definition: HDF5_File.h:500
Structure with 4D dimension sizes (3 in space and 1 in time).
hid_t OpenDataset(const hid_t parentGroup, TMatrixName &datasetName)
Open the HDF5 dataset at a specified place in the file tree.
Definition: HDF5_File.cpp:268
virtual ~TCuboidOutputHDF5Stream()
Destructor.