00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00032 #include "PetscTools.hpp"
00033 #include "Hdf5DataWriter.hpp"
00034
00035
00036 using std::string;
00037
00043 Hdf5DataWriter::Hdf5DataWriter(string directory, string baseName, bool cleanDirectory) :
00044 mDirectory(directory),
00045 mBaseName(baseName),
00046 mCleanDirectory(cleanDirectory),
00047 mIsInDefineMode(true),
00048 mIsFixedDimensionSet(false),
00049 mIsUnlimitedDimensionSet(false),
00050 mFileFixedDimensionSize(0U),
00051 mDataFixedDimensionSize(0U),
00052 mLo(0U),
00053 mHi(0U),
00054 mNumberOwned(0U),
00055 mOffset(0U),
00056 mIsDataComplete(true),
00057 mNeedExtend(false),
00058 mCurrentTimeStep(0)
00059 {
00060 int my_rank;
00061
00062 MPI_Comm_rank(PETSC_COMM_WORLD, &my_rank);
00063 if (my_rank==0)
00064 {
00065 mAmMaster=true;
00066 }
00067 else
00068 {
00069 mAmMaster=false;
00070 }
00071 }
00072
00073 Hdf5DataWriter::~Hdf5DataWriter()
00074 {
00075 }
00076
00084 void Hdf5DataWriter::DefineFixedDimension(long dimensionSize)
00085 {
00086 if (!mIsInDefineMode)
00087 {
00088 EXCEPTION("Cannot define variables when not in Define mode");
00089 }
00090 if (dimensionSize < 1)
00091 {
00092 EXCEPTION("Fixed dimension must be at least 1 long");
00093 }
00094 if (mIsFixedDimensionSet)
00095 {
00096 EXCEPTION("Fixed dimension already set");
00097 }
00098
00099
00100 mLo = DistributedVector::Begin().Global;
00101 mHi = DistributedVector::End().Global;
00102 mNumberOwned=mHi-mLo;
00103 mOffset=mLo;
00104 mFileFixedDimensionSize = dimensionSize;
00105 mDataFixedDimensionSize = dimensionSize;
00106 mIsFixedDimensionSet = true;
00107 }
00108
00117 void Hdf5DataWriter::DefineFixedDimension(std::vector<unsigned> nodesToOuput, long vecSize)
00118 {
00119 unsigned vector_size = nodesToOuput.size();
00120
00121 for (unsigned index=0; index < vector_size-1; index++)
00122 {
00123 if (nodesToOuput[index] >= nodesToOuput[index+1])
00124 {
00125 EXCEPTION("Input should be monotonic increasing");
00126 }
00127 }
00128
00129 if ((int) nodesToOuput.back() >= vecSize)
00130 {
00131 EXCEPTION("Vector size doesn't match nodes to ouput");
00132 }
00133
00134 DefineFixedDimension(vecSize);
00135
00136 mFileFixedDimensionSize = vector_size;
00137 mIsDataComplete = false;
00138 mIncompleteNodeIndices = nodesToOuput;
00139
00140 mOffset=0;
00141 mNumberOwned=0;
00142
00143 for (unsigned i=0;i<mIncompleteNodeIndices.size();i++)
00144 {
00145 if (mIncompleteNodeIndices[i] < mLo)
00146 {
00147 mOffset++;
00148 }
00149 else if(mIncompleteNodeIndices[i] < mHi)
00150 {
00151 mNumberOwned++;
00152 }
00153 }
00154 }
00155
00165 int Hdf5DataWriter::DefineVariable(string variableName, string variableUnits)
00166 {
00167 if (!mIsInDefineMode)
00168 {
00169 EXCEPTION("Cannot define variables when not in Define mode");
00170 }
00171
00172 CheckVariableName(variableName);
00173 CheckUnitsName(variableUnits);
00174
00175
00176 for (unsigned index=0; index<mVariables.size(); index++)
00177 {
00178 if (mVariables[index].mVariableName == variableName)
00179 {
00180 EXCEPTION("Variable name already exists");
00181 }
00182 }
00183
00184 DataWriterVariable new_variable;
00185 new_variable.mVariableName = variableName;
00186 new_variable.mVariableUnits = variableUnits;
00187 int variable_id;
00188
00189
00190 mVariables.push_back(new_variable);
00191
00192
00193 variable_id = mVariables.size()-1;
00194
00195 return variable_id;
00196 }
00197
00198
00199 void Hdf5DataWriter::CheckVariableName(std::string name)
00200 {
00201 if (name.length() == 0)
00202 {
00203 EXCEPTION("Variable name not allowed: may not be blank.");
00204 }
00205 CheckUnitsName(name);
00206 }
00207
00208 void Hdf5DataWriter::CheckUnitsName(std::string name)
00209 {
00210 for (unsigned i=0; i<name.length(); i++)
00211 {
00212 if (!isalnum(name[i]) && !(name[i]=='_'))
00213 {
00214 std::string error = "Variable name/units '" + name + "' not allowed: may only contain alphanumeric characters or '_'.";
00215 EXCEPTION(error);
00216 }
00217 }
00218 }
00219
00220 void Hdf5DataWriter::EndDefineMode()
00221 {
00222
00223 if (mVariables.size() < 1)
00224 {
00225 EXCEPTION("Cannot end define mode. No variables have been defined.");
00226 }
00227
00228
00229 if (mIsFixedDimensionSet == false)
00230 {
00231 EXCEPTION("Cannot end define mode. One fixed dimension should be defined.");
00232 }
00233
00234 mIsInDefineMode = false;
00235
00236 OutputFileHandler output_file_handler(mDirectory, mCleanDirectory);
00237 std::string results_dir = output_file_handler.GetOutputDirectoryFullPath();
00238 std::string file_name = results_dir + mBaseName + ".h5";
00239
00240
00241 hid_t property_list_id = H5Pcreate(H5P_FILE_ACCESS);
00242 H5Pset_fapl_mpio(property_list_id, PETSC_COMM_WORLD, MPI_INFO_NULL);
00243
00244
00245 mFileId = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, property_list_id);
00246 H5Pclose(property_list_id);
00247
00248
00249
00250
00251
00252
00253 mDatasetDims[0] = 1;
00254 mDatasetDims[1] = mFileFixedDimensionSize;
00255 mDatasetDims[2] = mVariables.size();
00256
00257 hsize_t* max_dims=NULL;
00258 hsize_t dataset_max_dims[DATASET_DIMS];
00259
00260 hid_t cparms=H5P_DEFAULT;
00261
00262 if (mIsUnlimitedDimensionSet)
00263 {
00264 dataset_max_dims[0] = H5S_UNLIMITED;
00265 dataset_max_dims[1] = mDatasetDims[1];
00266 dataset_max_dims[2] = mDatasetDims[2];
00267 max_dims = dataset_max_dims;
00268
00269
00270 hsize_t chunk_dims[DATASET_DIMS] ={1, mDatasetDims[1], mDatasetDims[2]};
00271 cparms = H5Pcreate (H5P_DATASET_CREATE);
00272 H5Pset_chunk( cparms, DATASET_DIMS, chunk_dims);
00273 }
00274
00275 hid_t filespace = H5Screate_simple(DATASET_DIMS, mDatasetDims, max_dims);
00276
00277
00278 mDatasetId = H5Dcreate(mFileId, "Data", H5T_NATIVE_DOUBLE, filespace, cparms);
00279 H5Sclose(filespace);
00280
00281
00282 const unsigned MAX_STRING_SIZE=100;
00283 hsize_t columns[1] = {mVariables.size()};
00284 hid_t colspace = H5Screate_simple(1, columns, NULL);
00285
00286
00287 char* col_data = (char*) malloc(mVariables.size() * sizeof(char) * MAX_STRING_SIZE);
00288
00289 char* col_data_offset = col_data;
00290 for (unsigned var=0; var<mVariables.size(); var++)
00291 {
00292 std::string full_name = mVariables[var].mVariableName + "("+mVariables[var].mVariableUnits+")";
00293 strcpy (col_data_offset, full_name.c_str());
00294 col_data_offset += sizeof(char) * MAX_STRING_SIZE;
00295 }
00296
00297
00298 hid_t string_type = H5Tcopy(H5T_C_S1);
00299 H5Tset_size(string_type, MAX_STRING_SIZE );
00300 hid_t attr = H5Acreate(mDatasetId, "Variable Details", string_type, colspace, H5P_DEFAULT );
00301
00302 H5Awrite(attr, string_type, col_data);
00303
00304
00305 free(col_data);
00306 H5Sclose(colspace);
00307 H5Aclose(attr);
00308
00309
00310 columns[0] = 1;
00311 colspace = H5Screate_simple(1, columns, NULL);
00312 attr = H5Acreate(mDatasetId, "IsDataComplete", H5T_NATIVE_UINT, colspace, H5P_DEFAULT );
00313
00314
00315 unsigned is_data_complete= mIsDataComplete ? 1 : 0;
00316 H5Awrite(attr, H5T_NATIVE_UINT, &is_data_complete);
00317
00318 H5Sclose(colspace);
00319 H5Aclose(attr);
00320
00321 if (!mIsDataComplete)
00322 {
00323
00324
00325 columns[0] = mFileFixedDimensionSize;
00326 colspace = H5Screate_simple(1, columns, NULL);
00327 attr = H5Acreate(mDatasetId, "NodeMap", H5T_NATIVE_UINT, colspace, H5P_DEFAULT );
00328
00329
00330 H5Awrite(attr, H5T_NATIVE_UINT, &mIncompleteNodeIndices[0]);
00331
00332 H5Sclose(colspace);
00333 H5Aclose(attr);
00334 }
00335
00336
00337
00338
00339
00340 if (mIsUnlimitedDimensionSet)
00341 {
00342 hsize_t time_dataset_dims[1] = {1};
00343 hsize_t time_dataset_max_dims[1] = {H5S_UNLIMITED};
00344
00345
00346 hsize_t time_chunk_dims[1] ={1};
00347 hid_t time_cparms = H5Pcreate (H5P_DATASET_CREATE);
00348 H5Pset_chunk( time_cparms, 1, time_chunk_dims);
00349
00350 hid_t time_filespace = H5Screate_simple(1, time_dataset_dims, time_dataset_max_dims);
00351
00352
00353 mTimeDatasetId = H5Dcreate(mFileId, mUnlimitedDimensionName.c_str(), H5T_NATIVE_DOUBLE, time_filespace, time_cparms);
00354
00355
00356 hsize_t one = 1;
00357 hid_t one_column_space = H5Screate_simple(1, &one, NULL);
00358
00359
00360 hid_t unit_attr = H5Acreate(mTimeDatasetId, "Unit", string_type, one_column_space, H5P_DEFAULT);
00361
00362
00363 char unit_string[MAX_STRING_SIZE];
00364 strcpy(unit_string, mUnlimitedDimensionUnit.c_str());
00365 H5Awrite(unit_attr, string_type, unit_string);
00366
00367
00368 H5Sclose(one_column_space);
00369 H5Aclose(unit_attr);
00370 H5Sclose(time_filespace);
00371 }
00372
00373 }
00374
00375 void Hdf5DataWriter::PutVector(int variableID, Vec petscVector)
00376 {
00377 if (mIsInDefineMode)
00378 {
00379 EXCEPTION("Cannot write data while in define mode.");
00380 }
00381
00382 int vector_size;
00383 VecGetSize(petscVector, &vector_size);
00384
00385 if ((unsigned) vector_size != mDataFixedDimensionSize)
00386 {
00387 EXCEPTION("Vector size doesn't match fixed dimension");
00388 }
00389
00390
00391 PossiblyExtend();
00392
00393
00394 hid_t memspace=0;
00395 if (mNumberOwned !=0)
00396 {
00397 hsize_t v_size[1] = {mNumberOwned};
00398 memspace = H5Screate_simple(1, v_size, NULL);
00399 }
00400
00401 hsize_t count[DATASET_DIMS] = {1, mNumberOwned, 1};
00402 hsize_t offset_dims[DATASET_DIMS] = {mCurrentTimeStep, mOffset, variableID};
00403 hid_t file_dataspace = H5Dget_space(mDatasetId);
00404
00405
00406 hid_t property_list_id = H5Pcreate(H5P_DATASET_XFER);
00407 H5Pset_dxpl_mpio(property_list_id, H5FD_MPIO_COLLECTIVE);
00408
00409 H5Sselect_hyperslab(file_dataspace, H5S_SELECT_SET, offset_dims, NULL, count, NULL);
00410
00411 double* p_petsc_vector;
00412 VecGetArray(petscVector, &p_petsc_vector);
00413
00414 if (mIsDataComplete)
00415 {
00416 H5Dwrite(mDatasetId, H5T_NATIVE_DOUBLE, memspace, file_dataspace, property_list_id, p_petsc_vector);
00417 }
00418 else
00419 {
00420
00421 double local_data[mNumberOwned];
00422 for (unsigned i=0;i<mNumberOwned;i++)
00423 {
00424 local_data[i] = p_petsc_vector[ mIncompleteNodeIndices[mOffset+i]-mLo ];
00425
00426 }
00427 H5Dwrite(mDatasetId, H5T_NATIVE_DOUBLE, memspace, file_dataspace, property_list_id, local_data);
00428 }
00429
00430 VecRestoreArray(petscVector, &p_petsc_vector);
00431
00432 H5Pclose(property_list_id);
00433 H5Sclose(file_dataspace);
00434 if (mNumberOwned !=0)
00435 {
00436 H5Sclose(memspace);
00437 }
00438 }
00439
00440 void Hdf5DataWriter::PutStripedVector(int firstVariableID, int secondVariableID, Vec petscVector)
00441 {
00442 if (mIsInDefineMode)
00443 {
00444 EXCEPTION("Cannot write data while in define mode.");
00445 }
00446
00447 int NUM_STRIPES=2;
00448
00449
00450 if (secondVariableID-firstVariableID != 1)
00451 {
00452 EXCEPTION("Columns should be consecutive. Try reordering them.");
00453 }
00454
00455 int vector_size;
00456 VecGetSize(petscVector, &vector_size);
00457
00458 if ((unsigned) vector_size != NUM_STRIPES*mDataFixedDimensionSize)
00459 {
00460 EXCEPTION("Vector size doesn't match fixed dimension");
00461 }
00462
00463
00464 PossiblyExtend();
00465
00466
00467 hid_t memspace=0;
00468 if (mNumberOwned !=0)
00469 {
00470 hsize_t v_size[1] = {mNumberOwned*NUM_STRIPES};
00471 memspace = H5Screate_simple(1, v_size, NULL);
00472 }
00473
00474
00475 hsize_t start[DATASET_DIMS] = {mCurrentTimeStep, mOffset, firstVariableID};
00476 hsize_t stride[DATASET_DIMS] = {1, 1, secondVariableID-firstVariableID};
00477 hsize_t block_size[DATASET_DIMS] = {1, mNumberOwned, 1};
00478 hsize_t number_blocks[DATASET_DIMS] = {1, 1, NUM_STRIPES};
00479
00480 hid_t hyperslab_space = H5Dget_space(mDatasetId);
00481 H5Sselect_hyperslab(hyperslab_space, H5S_SELECT_SET, start, stride, number_blocks, block_size);
00482
00483
00484 hid_t property_list_id = H5Pcreate(H5P_DATASET_XFER);
00485 H5Pset_dxpl_mpio(property_list_id, H5FD_MPIO_COLLECTIVE);
00486
00487 double* p_petsc_vector;
00488 VecGetArray(petscVector, &p_petsc_vector);
00489
00490 if (mIsDataComplete)
00491 {
00492 H5Dwrite(mDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, p_petsc_vector);
00493 }
00494 else
00495 {
00496
00497 double local_data[mNumberOwned*NUM_STRIPES];
00498 for (unsigned i=0;i<mNumberOwned;i++)
00499 {
00501 unsigned local_node_number=mIncompleteNodeIndices[mOffset+i] - mLo;
00502 local_data[NUM_STRIPES*i] = p_petsc_vector[ local_node_number*NUM_STRIPES ];
00503 local_data[NUM_STRIPES*i+1] = p_petsc_vector[ local_node_number*NUM_STRIPES + 1];
00504 }
00505 H5Dwrite(mDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, local_data);
00506 }
00507
00508 VecRestoreArray(petscVector, &p_petsc_vector);
00509
00510 H5Sclose(hyperslab_space);
00511 if (mNumberOwned !=0)
00512 {
00513 H5Sclose(memspace);
00514 }
00515 H5Pclose(property_list_id);
00516 }
00517
00518 void Hdf5DataWriter::PutUnlimitedVariable(double value)
00519 {
00520 if (mIsInDefineMode)
00521 {
00522 EXCEPTION("Cannot write data while in define mode.");
00523 }
00524
00525 if(!mIsUnlimitedDimensionSet)
00526 {
00527 EXCEPTION("PutUnlimitedVariable() called but no unlimited dimension has been set");
00528 }
00529
00530
00531 if(!mAmMaster)
00532 {
00533 return;
00534 }
00535
00536
00537 PossiblyExtend();
00538
00539 hsize_t size[1] = {1};
00540 hid_t memspace = H5Screate_simple(1, size, NULL);
00541
00542
00543 hsize_t count[1] = {1};
00544 hsize_t offset[1] = {mCurrentTimeStep};
00545 hid_t hyperslab_space = H5Dget_space(mTimeDatasetId);
00546 H5Sselect_hyperslab(hyperslab_space, H5S_SELECT_SET, offset, NULL, count, NULL);
00547
00548 H5Dwrite(mTimeDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, H5P_DEFAULT, &value);
00549
00550 H5Sclose(hyperslab_space);
00551 H5Sclose(memspace);
00552 }
00553
00554
00555 void Hdf5DataWriter::Close()
00556 {
00557 if (mIsInDefineMode)
00558 {
00559 return;
00560 }
00561 H5Dclose(mDatasetId);
00562
00563 if (mIsUnlimitedDimensionSet)
00564 {
00565 H5Dclose(mTimeDatasetId);
00566 }
00567
00568 H5Fclose(mFileId);
00569 }
00570
00571
00572 void Hdf5DataWriter::DefineUnlimitedDimension(std::string variableName, std::string variableUnits)
00573 {
00574 if (mIsUnlimitedDimensionSet)
00575 {
00576 EXCEPTION("Unlimited dimension already set. Cannot be defined twice");
00577 }
00578
00579 if (!mIsInDefineMode)
00580 {
00581 EXCEPTION("Cannot define variables when not in Define mode");
00582 }
00583
00584 mIsUnlimitedDimensionSet = true;
00585 mUnlimitedDimensionName = variableName;
00586 mUnlimitedDimensionUnit = variableUnits;
00587 }
00588
00589 void Hdf5DataWriter::AdvanceAlongUnlimitedDimension()
00590 {
00591 if (!mIsUnlimitedDimensionSet)
00592 {
00593 EXCEPTION("Trying to advance along an unlimited dimension without having defined any");
00594 }
00595
00596
00597 mDatasetDims[0]++;
00598 mNeedExtend=true;
00599
00600 mCurrentTimeStep++;
00601 }
00602
00603 void Hdf5DataWriter::PossiblyExtend()
00604 {
00605 if (mNeedExtend)
00606 {
00607 H5Dextend (mDatasetId, mDatasetDims);
00608 H5Dextend (mTimeDatasetId, mDatasetDims);
00609 }
00610 mNeedExtend=false;
00611 }