XdmfHDF5Writer.cpp 16.3 KB
Newer Older
Kenneth Leiter's avatar
Kenneth Leiter committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
/*****************************************************************************/
/*                                    XDMF                                   */
/*                       eXtensible Data Model and Format                    */
/*                                                                           */
/*  Id : XdmfHDF5Writer.cpp                                                  */
/*                                                                           */
/*  Author:                                                                  */
/*     Kenneth Leiter                                                        */
/*     kenneth.leiter@arl.army.mil                                           */
/*     US Army Research Laboratory                                           */
/*     Aberdeen Proving Ground, MD                                           */
/*                                                                           */
/*     Copyright @ 2011 US Army Research Laboratory                          */
/*     All Rights Reserved                                                   */
/*     See Copyright.txt for details                                         */
/*                                                                           */
/*     This software is distributed WITHOUT ANY WARRANTY; without            */
/*     even the implied warranty of MERCHANTABILITY or FITNESS               */
/*     FOR A PARTICULAR PURPOSE.  See the above copyright notice             */
/*     for more information.                                                 */
/*                                                                           */
/*****************************************************************************/
23

24
#include <hdf5.h>
25
#include <cmath>
26
#include <cstdio>
27 28
#include <numeric>
#include <sstream>
29
#include "XdmfArray.hpp"
30
#include "XdmfArrayType.hpp"
31
#include "XdmfError.hpp"
32
#include "XdmfHDF5Controller.hpp"
33
#include "XdmfHDF5Writer.hpp"
34

35 36
namespace {

37
  const static unsigned int DEFAULT_CHUNK_SIZE = 1000;
38 39 40

}

41 42 43 44 45 46 47 48
/**
 * PIMPL
 */
class XdmfHDF5Writer::XdmfHDF5WriterImpl {

public:

  XdmfHDF5WriterImpl():
49
    mChunkSize(DEFAULT_CHUNK_SIZE),
50 51 52 53 54 55 56 57 58 59 60 61 62
    mHDF5Handle(-1)
  {
  };

  ~XdmfHDF5WriterImpl()
  {
    closeFile();
  };

  void
  closeFile()
  {
    if(mHDF5Handle >= 0) {
63
      /*herr_t status =*/H5Fclose(mHDF5Handle);
64 65
      mHDF5Handle = -1;
    }
66 67
  };
  
68 69
  int
  openFile(const std::string & filePath,
70 71
           const int fapl,
           const int mDataSetId)
72 73 74 75 76 77 78 79 80 81 82 83
  {
    if(mHDF5Handle >= 0) {
      // Perhaps we should throw a warning.
      closeFile();
    }

    // Save old error handler and turn off error handling for now
    H5E_auto_t old_func;
    void * old_client_data;
    H5Eget_auto(0, &old_func, &old_client_data);
    H5Eset_auto2(0, NULL, NULL);
  
84
    int toReturn = 0;
85 86 87 88 89

    if(H5Fis_hdf5(filePath.c_str()) > 0) {
      mHDF5Handle = H5Fopen(filePath.c_str(), 
                            H5F_ACC_RDWR, 
                            fapl);
90 91 92 93 94 95 96 97 98
      if(mDataSetId == 0) {
        hsize_t numObjects;
        /*herr_t status = */H5Gget_num_objs(mHDF5Handle,
                                            &numObjects);
        toReturn = numObjects;
      }
      else {
        toReturn = mDataSetId;
      }
99 100 101 102 103 104 105 106 107 108 109
    }
    else {
      mHDF5Handle = H5Fcreate(filePath.c_str(),
                              H5F_ACC_TRUNC,
                              H5P_DEFAULT,
                              fapl);
    }

    // Restore previous error handler
    H5Eset_auto2(0, old_func, old_client_data);

110
    return toReturn;
111 112
  }

113
  unsigned int mChunkSize;
114 115 116
  hid_t mHDF5Handle;

};
117

118
shared_ptr<XdmfHDF5Writer>
119 120
XdmfHDF5Writer::New(const std::string & filePath,
                    const bool clobberFile)
121
{
122 123 124
  if(clobberFile) {
    std::remove(filePath.c_str());
  }
125
  shared_ptr<XdmfHDF5Writer> p(new XdmfHDF5Writer(filePath));
126
  return p;
127 128
}

129
XdmfHDF5Writer::XdmfHDF5Writer(const std::string & filePath) :
130 131
  XdmfHeavyDataWriter(filePath),
  mImpl(new XdmfHDF5WriterImpl())
132 133 134 135 136
{
}

XdmfHDF5Writer::~XdmfHDF5Writer()
{
137
  delete mImpl;
138 139
}

140
shared_ptr<XdmfHDF5Controller>
141 142
XdmfHDF5Writer::createHDF5Controller(const std::string & hdf5FilePath,
                                     const std::string & dataSetPath,
143
                                     const shared_ptr<const XdmfArrayType> type,
144 145
                                     const std::vector<unsigned int> & start,
                                     const std::vector<unsigned int> & stride,
146 147
                                     const std::vector<unsigned int> & dimensions,
                                     const std::vector<unsigned int> & dataspaceDimensions)
148
{
149 150 151 152 153
  return XdmfHDF5Controller::New(hdf5FilePath,
                                 dataSetPath,
                                 type,
                                 start,
                                 stride,
154 155
                                 dimensions,
                                 dataspaceDimensions);
156 157
}

158 159 160 161 162 163
void 
XdmfHDF5Writer::closeFile()
{
  mImpl->closeFile();
}

164 165 166 167 168 169
unsigned int
XdmfHDF5Writer::getChunkSize() const
{
  return mImpl->mChunkSize;
}

170 171 172
void 
XdmfHDF5Writer::openFile()
{
173 174 175 176 177 178 179
  this->openFile(H5P_DEFAULT);
}

void
XdmfHDF5Writer::openFile(const int fapl)
{
  mDataSetId = mImpl->openFile(mFilePath,
180 181
                               fapl,
                               mDataSetId);
182 183
}

184 185 186 187 188 189
void
XdmfHDF5Writer::setChunkSize(const unsigned int chunkSize)
{
  mImpl->mChunkSize = chunkSize;
}

190 191
void
XdmfHDF5Writer::visit(XdmfArray & array,
192
                      const shared_ptr<XdmfBaseVisitor> visitor)
193
{
194
  this->write(array, H5P_DEFAULT);
195 196
}

197 198 199
void
XdmfHDF5Writer::write(XdmfArray & array,
                      const int fapl)
200
{
201
  hid_t datatype = -1;
Kenneth Leiter's avatar
Kenneth Leiter committed
202
  bool closeDatatype = false;
203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231

  if(array.isInitialized()) {
    if(array.getArrayType() == XdmfArrayType::Int8()) {
      datatype = H5T_NATIVE_CHAR;
    }
    else if(array.getArrayType() == XdmfArrayType::Int16()) {
      datatype = H5T_NATIVE_SHORT;
    }
    else if(array.getArrayType() == XdmfArrayType::Int32()) {
      datatype = H5T_NATIVE_INT;
    }
    else if(array.getArrayType() == XdmfArrayType::Int64()) {
      datatype = H5T_NATIVE_LONG;
    }
    else if(array.getArrayType() == XdmfArrayType::Float32()) {
      datatype = H5T_NATIVE_FLOAT;
    }
    else if(array.getArrayType() == XdmfArrayType::Float64()) {
      datatype = H5T_NATIVE_DOUBLE;
    }
    else if(array.getArrayType() == XdmfArrayType::UInt8()) {
      datatype = H5T_NATIVE_UCHAR;
    }
    else if(array.getArrayType() == XdmfArrayType::UInt16()) {
      datatype = H5T_NATIVE_USHORT;
    }
    else if(array.getArrayType() == XdmfArrayType::UInt32()) {
      datatype = H5T_NATIVE_UINT;
    }
Kenneth Leiter's avatar
Kenneth Leiter committed
232 233 234 235 236
    else if(array.getArrayType() == XdmfArrayType::String()) {
      datatype = H5Tcopy(H5T_C_S1);
      H5Tset_size(datatype, H5T_VARIABLE);
      closeDatatype = true;
    }
237
    else {
238 239 240
      XdmfError::message(XdmfError::FATAL,
                         "Array of unsupported type in "
                         "XdmfHDF5Writer::write");
241 242 243 244 245
    }
  }

  if(datatype != -1) {
    std::string hdf5FilePath = mFilePath;
246

247 248
    std::stringstream dataSetPath;

249 250 251 252 253 254 255 256 257 258
    shared_ptr<XdmfHeavyDataController> heavyDataController = 
      array.getHeavyDataController();
    const std::vector<unsigned int> & dimensions = array.getDimensions();
    std::vector<unsigned int> dataspaceDimensions = dimensions;
    std::vector<unsigned int> start(dimensions.size(), 0);
    std::vector<unsigned int> stride(dimensions.size(), 1);

    if((mMode == Overwrite || mMode == Append || mMode == Hyperslab)
       && heavyDataController) {
      
259
      // Write to the previous dataset
260 261 262 263 264 265 266
      dataSetPath << heavyDataController->getDataSetPath();
      hdf5FilePath = heavyDataController->getFilePath();
      if(mMode == Hyperslab) {
        dataspaceDimensions = heavyDataController->getDataspaceDimensions();
        start = heavyDataController->getStart();
        stride = heavyDataController->getStride();
      }
267 268 269 270 271 272 273 274 275 276 277 278 279
    }
    else {
      dataSetPath << "Data" << mDataSetId;
    }

    // Open a hdf5 dataset and write to it on disk.
    herr_t status;

    // Save old error handler and turn off error handling for now
    H5E_auto_t old_func;
    void * old_client_data;
    H5Eget_auto(0, &old_func, &old_client_data);
    H5Eset_auto2(0, NULL, NULL);
280 281 282
   
    bool closeFile = false;
    if(mImpl->mHDF5Handle < 0) {
283
      mImpl->openFile(hdf5FilePath,
284 285
                      fapl,
                      mDataSetId);
286
      closeFile = true;
287
    }
288 289

    hid_t dataset = H5Dopen(mImpl->mHDF5Handle,
290 291
                            dataSetPath.str().c_str(),
                            H5P_DEFAULT);
292 293
    // if default mode find a new data set to write to (keep
    // incrementing dataSetId)
294
    if(dataset >= 0 && mMode == Default) {
295
      H5Dclose(dataset);
296 297 298 299 300 301 302 303 304 305 306 307
      while(true) {
        dataSetPath.str(std::string());
        dataSetPath << "Data" << ++mDataSetId;
        if(!H5Lexists(mImpl->mHDF5Handle,
                      dataSetPath.str().c_str(),
                      H5P_DEFAULT)) {
          dataset = H5Dopen(mImpl->mHDF5Handle,
                            dataSetPath.str().c_str(),
                            H5P_DEFAULT); 
          break;
        }
      }
308
    }
309

310 311 312
    // Restore previous error handler
    H5Eset_auto2(0, old_func, old_client_data);

313 314 315
    hid_t dataspace = H5S_ALL;
    hid_t memspace = H5S_ALL;

316 317
    std::vector<hsize_t> current_dims(dataspaceDimensions.begin(), 
                                      dataspaceDimensions.end());
318

319
    if(dataset < 0) {
320
      std::vector<hsize_t> maximum_dims(dimensions.size(), H5S_UNLIMITED);
321 322 323
      dataspace = H5Screate_simple(dimensions.size(),
                                   &current_dims[0],
                                   &maximum_dims[0]);
324 325 326

      // calculate a proper chunk size - for multidimensional datasets the
      // chunk dimensions have similar shape to original dataset
327
      hid_t property = H5Pcreate(H5P_DATASET_CREATE);
328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344
      const hsize_t totalDimensionsSize = 
        std::accumulate(current_dims.begin(),
                        current_dims.end(),
                        1,
                        std::multiplies<hsize_t>());
      const double factor = 
        std::pow(((double)mImpl->mChunkSize / totalDimensionsSize), 
                 1.0 / current_dims.size());
      std::vector<hsize_t> chunk_size(current_dims.begin(),
                                      current_dims.end());
      for(std::vector<hsize_t>::iterator iter = chunk_size.begin();
          iter != chunk_size.end(); ++iter) {
        *iter = (hsize_t)(*iter * factor);
        if(*iter == 0) {
          *iter = 1;
        }
      }
345
      status = H5Pset_chunk(property, dimensions.size(), &chunk_size[0]);
346
      dataset = H5Dcreate(mImpl->mHDF5Handle,
347 348
                          dataSetPath.str().c_str(),
                          datatype,
349
                          dataspace,
350 351 352
                          H5P_DEFAULT,
                          property,
                          H5P_DEFAULT);
353
      if(dataset < 0) {
354 355 356
        std::stringstream error;
        error << "H5Dcreate returned failure in XdmfHDF5Writer::write -- " 
              << "dataset: " << dataset;
357
        XdmfError::message(XdmfError::FATAL,
358
                           error.str());
359
      }
360 361
      status = H5Pclose(property);
    }
362 363

    if(mMode == Append) {
364
      // Need to resize dataset to fit new data
365 366 367 368 369 370 371
      
      // Get size of old dataset
      dataspace = H5Dget_space(dataset);
      hssize_t datasize = H5Sget_simple_extent_npoints(dataspace);
      status = H5Sclose(dataspace);
      
      // Resize to fit size of old and new data.
372
      hsize_t size = array.getSize();
373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432
      hsize_t newSize = size + datasize;
      status = H5Dset_extent(dataset, &newSize);
      
      // Select hyperslab to write to.
      memspace = H5Screate_simple(1, &size, NULL);
      dataspace = H5Dget_space(dataset); 
      hsize_t start = datasize;
      status = H5Sselect_hyperslab(dataspace,
                                   H5S_SELECT_SET,
                                   &start,
                                   NULL,
                                   &size,
                                   NULL);
    }
    else if(mMode == Overwrite) {
      // Overwriting - dataset rank must remain the same (hdf5 constraint)
      dataspace = H5Dget_space(dataset);
      
      const unsigned int ndims = H5Sget_simple_extent_ndims(dataspace);
      if(ndims != current_dims.size())
        XdmfError::message(XdmfError::FATAL,                            \
                           "Data set rank different -- ndims != "
                           "current_dims.size() -- in "
                           "XdmfHDF5Writer::write");

      status = H5Dset_extent(dataset, &current_dims[0]);
      dataspace = H5Dget_space(dataset);
    }
    else if(mMode == Hyperslab) {
      // Hyperslab - dataset rank must remain the same (hdf5 constraint)
      dataspace = H5Dget_space(dataset);
      
      const unsigned int ndims = H5Sget_simple_extent_ndims(dataspace);
      if(ndims != current_dims.size())
        XdmfError::message(XdmfError::FATAL,                            \
                           "Data set rank different -- ndims != "
                           "current_dims.size() -- in "
                           "XdmfHDF5Writer::write");
      
      status = H5Dset_extent(dataset, &current_dims[0]);
      dataspace = H5Dget_space(dataset);

      std::vector<hsize_t> count(dimensions.begin(),
                                 dimensions.end());
      std::vector<hsize_t> currStride(stride.begin(),
                                      stride.end());
      std::vector<hsize_t> currStart(start.begin(),
                                     start.end());

      memspace = H5Screate_simple(count.size(),
                                  &(count[0]),
                                  NULL);
      status = H5Sselect_hyperslab(dataspace,
                                   H5S_SELECT_SET,
                                   &currStart[0],
                                   &currStride[0],
                                   &count[0],
                                   NULL) ;
      
      if(status < 0) {
433 434 435
        std::stringstream error;
        error << "H5Dset_extent returned failure in XdmfHDF5Writer::write "
              << "-- status: " << status;
436
        XdmfError::message(XdmfError::FATAL,
437
                           error.str());
438 439
      }
    }
440

441 442 443 444 445 446
    status = H5Dwrite(dataset,
                      datatype,
                      memspace,
                      dataspace,
                      H5P_DEFAULT,
                      array.getValuesInternal());
447 448

    if(status < 0) {
449 450 451
      std::stringstream error;
      error << "H5Dwrite returned failure in XdmfHDF5Writer::write "
            << "-- status: " << status;
452
      XdmfError::message(XdmfError::FATAL, 
453
                         error.str());
454 455
    }

456 457 458 459 460 461 462
    if(dataspace != H5S_ALL) {
      status = H5Sclose(dataspace);
    }
    if(memspace != H5S_ALL) {
      status = H5Sclose(memspace);
    }
    status = H5Dclose(dataset);
Kenneth Leiter's avatar
Kenneth Leiter committed
463 464 465
    if(closeDatatype) {
      status = H5Tclose(datatype);
    }
466 467 468
    if(closeFile) {
      mImpl->closeFile();
    }
469

470
    // Attach a new controller to the array
471 472
    shared_ptr<XdmfHDF5Controller> newDataController =
      shared_ptr<XdmfHDF5Controller>();
473

474
    unsigned int newSize = array.getSize();
475 476
    if(mMode == Append && heavyDataController) {
      newSize = array.getSize() + heavyDataController->getSize();
477 478 479 480 481 482
      newDataController =
        this->createHDF5Controller(hdf5FilePath,
                                   dataSetPath.str(),
                                   array.getArrayType(),
                                   std::vector<unsigned int>(1, 0),
                                   std::vector<unsigned int>(1, 1),
483
                                   std::vector<unsigned int>(1, newSize),
484 485
                                   std::vector<unsigned int>(1, newSize));
    }
486

487
    if(mMode == Default || !heavyDataController) {
488
      ++mDataSetId;
489
    }
490

491 492 493 494 495
    if(!newDataController) {
      newDataController =
        this->createHDF5Controller(hdf5FilePath,
                                   dataSetPath.str(),
                                   array.getArrayType(),
496 497 498 499
                                   start,
                                   stride,
                                   dimensions,
                                   dataspaceDimensions);
500
    }
501
    array.setHeavyDataController(newDataController);
502 503 504 505

    if(mReleaseData) {
      array.release();
    }
506
  }
507
}