Hdf5DataWriter.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 /*
00029 * Implementation file for Hdf5DataWriter class.
00030 *
00031 */
00032 #include <set>
00033 #include <cstring> //For strcmp etc. Needed in gcc-4.4
00034 
00035 #include "Hdf5DataWriter.hpp"
00036 
00037 #include "Exception.hpp"
00038 #include "OutputFileHandler.hpp"
00039 #include "PetscTools.hpp"
00040 #include "Version.hpp"
00041 
00042 Hdf5DataWriter::Hdf5DataWriter(DistributedVectorFactory& rVectorFactory,
00043                                const std::string& rDirectory,
00044                                const std::string& rBaseName,
00045                                bool cleanDirectory,
00046                                bool extendData)
00047     : mrVectorFactory(rVectorFactory),
00048       mDirectory(rDirectory),
00049       mBaseName(rBaseName),
00050       mCleanDirectory(cleanDirectory),
00051       mIsInDefineMode(true),
00052       mIsFixedDimensionSet(false),
00053       mIsUnlimitedDimensionSet(false),
00054       mEstimatedUnlimitedLength(1u),
00055       mFileFixedDimensionSize(0u),
00056       mDataFixedDimensionSize(0u),
00057       mLo(mrVectorFactory.GetLow()),
00058       mHi(mrVectorFactory.GetHigh()),
00059       mNumberOwned(0u),
00060       mOffset(0u),
00061       mIsDataComplete(true),
00062       mNeedExtend(false),
00063       mUseMatrixForIncompleteData(false),
00064       mCurrentTimeStep(0),
00065       mSinglePermutation(NULL),
00066       mDoublePermutation(NULL),
00067       mSingleIncompleteOutputMatrix(NULL),
00068       mDoubleIncompleteOutputMatrix(NULL),
00069       mUseOptimalChunkSizeAlgorithm(true),
00070       mFixedChunkSize(0)
00071 {
00072     if (extendData && cleanDirectory)
00073     {
00074         EXCEPTION("You are asking to delete a file and then extend it, change arguments to constructor.");
00075     }
00076 
00077     if (extendData)
00078     {
00079         // Where to find the file
00080         OutputFileHandler output_file_handler(mDirectory, false);
00081         std::string results_dir = output_file_handler.GetOutputDirectoryFullPath();
00082         std::string file_name = results_dir + mBaseName + ".h5";
00083 
00084         // Set up a property list saying how we'll open the file
00085         hid_t property_list_id = H5Pcreate(H5P_FILE_ACCESS);
00086         H5Pset_fapl_mpio(property_list_id, PETSC_COMM_WORLD, MPI_INFO_NULL);
00087 
00088         // Open the file and free the property list
00089         mFileId = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, property_list_id);
00090         H5Pclose(property_list_id);
00091 
00092         if (mFileId < 0)
00093         {
00094             mDatasetId = 0;
00095             EXCEPTION("Hdf5DataWriter could not open " + file_name);
00096         }
00097 
00098         // Open the main dataset, and figure out its size/shape
00099         mDatasetId = H5Dopen(mFileId, "Data");
00100         hid_t variables_dataspace = H5Dget_space(mDatasetId);
00101         //unsigned variables_dataset_rank = H5Sget_simple_extent_ndims(variables_dataspace);
00102         hsize_t dataset_max_sizes[DATASET_DIMS];
00103         H5Sget_simple_extent_dims(variables_dataspace, mDatasetDims, dataset_max_sizes);
00104         H5Sclose(variables_dataspace);
00105         // Check that an unlimited dimension is defined
00106         if (dataset_max_sizes[0] != H5S_UNLIMITED)
00107         {
00108             H5Dclose(mDatasetId);
00109             H5Fclose(mFileId);
00110             EXCEPTION("Tried to open a datafile for extending which doesn't have an unlimited dimension.");
00111         }
00112         mIsUnlimitedDimensionSet = true;
00113         // Sanity check other dimension sizes
00114         for (unsigned i=1; i<DATASET_DIMS; i++)  // Zero is excluded since it is unlimited
00115         {
00116             assert(mDatasetDims[i] == dataset_max_sizes[i]);
00117         }
00118         mFileFixedDimensionSize = mDatasetDims[1];
00119         mIsFixedDimensionSet = true;
00120         mVariables.reserve(mDatasetDims[2]);
00121 
00122         // Figure out what the variables are
00123         hid_t attribute_id = H5Aopen_name(mDatasetId, "Variable Details");
00124         hid_t attribute_type  = H5Aget_type(attribute_id);
00125         hid_t attribute_space = H5Aget_space(attribute_id);
00126         hsize_t attr_dataspace_dim;
00127         H5Sget_simple_extent_dims(attribute_space, &attr_dataspace_dim, NULL);
00128         unsigned num_columns = H5Sget_simple_extent_npoints(attribute_space);
00129         assert(num_columns == mDatasetDims[2]); // I think...
00130 
00131         char* string_array = (char *)malloc(sizeof(char)*MAX_STRING_SIZE*(int)num_columns);
00132         H5Aread(attribute_id, attribute_type, string_array);
00133         // Loop over columns/variables
00134         for (unsigned index=0; index<num_columns; index++)
00135         {
00136             // Read the string from the raw vector
00137             std::string column_name_unit(&string_array[MAX_STRING_SIZE*index]);
00138             // Find location of unit name
00139             size_t name_length = column_name_unit.find('(');
00140             size_t unit_length = column_name_unit.find(')') - name_length - 1;
00141             // Create the variable
00142             DataWriterVariable var;
00143             var.mVariableName = column_name_unit.substr(0, name_length);
00144             var.mVariableUnits = column_name_unit.substr(name_length+1, unit_length);
00145             mVariables.push_back(var);
00146         }
00147         // Free memory, release ids
00148         free(string_array);
00149         H5Tclose(attribute_type);
00150         H5Sclose(attribute_space);
00151         H5Aclose(attribute_id);
00152 
00153         // Now deal with time
00154         mUnlimitedDimensionName = "Time"; // Assumed by the reader...
00155         mTimeDatasetId = H5Dopen(mFileId, mUnlimitedDimensionName.c_str());
00156         mUnlimitedDimensionUnit = "ms"; // Assumed by Chaste...
00157         // How many time steps have been written so far?
00158         hid_t timestep_dataspace = H5Dget_space(mTimeDatasetId);
00159         hsize_t num_timesteps;
00160         H5Sget_simple_extent_dims(timestep_dataspace, &num_timesteps, NULL);
00161         H5Sclose(timestep_dataspace);
00162         mCurrentTimeStep = (long)num_timesteps - 1;
00163 
00164         // Incomplete data?
00165         attribute_id = H5Aopen_name(mDatasetId, "IsDataComplete");
00166         if (attribute_id < 0)
00167         {
00168 #define COVERAGE_IGNORE
00169             // Old format, before we added the attribute.
00170             EXCEPTION("Extending old-format files isn't supported.");
00171 #undef COVERAGE_IGNORE
00172         }
00173         else
00174         {
00175             attribute_type = H5Aget_type(attribute_id);
00176             attribute_space = H5Aget_space(attribute_id);
00177             unsigned is_data_complete;
00178             H5Aread(attribute_id, H5T_NATIVE_UINT, &is_data_complete);
00179             mIsDataComplete = (is_data_complete == 1) ? true : false;
00180             H5Tclose(attribute_type);
00181             H5Sclose(attribute_space);
00182             H5Aclose(attribute_id);
00183         }
00184         if (mIsDataComplete)
00185         {
00186             mNumberOwned = mrVectorFactory.GetLocalOwnership();
00187             mOffset = mLo;
00188             mDataFixedDimensionSize = mFileFixedDimensionSize;
00189         }
00190         else
00191         {
00192             // Read which nodes appear in the file (mIncompleteNodeIndices)
00193             attribute_id = H5Aopen_name(mDatasetId, "NodeMap");
00194             attribute_type  = H5Aget_type(attribute_id);
00195             attribute_space = H5Aget_space(attribute_id);
00196             // Get the dataset/dataspace dimensions
00197             unsigned num_node_indices = H5Sget_simple_extent_npoints(attribute_space);
00198             // Read data from hyperslab in the file into the hyperslab in memory
00199             mIncompleteNodeIndices.clear();
00200             mIncompleteNodeIndices.resize(num_node_indices);
00201             H5Aread(attribute_id, H5T_NATIVE_UINT, &mIncompleteNodeIndices[0]);
00202             // Release ids
00203             H5Tclose(attribute_type);
00204             H5Sclose(attribute_space);
00205             H5Aclose(attribute_id);
00206             // Set up what data we can
00207             mNumberOwned = mrVectorFactory.GetLocalOwnership();
00208             ComputeIncompleteOffset();
00212             mDataFixedDimensionSize = UINT_MAX;
00213             H5Dclose(mDatasetId);
00214             H5Dclose(mTimeDatasetId);
00215             H5Fclose(mFileId);
00216             EXCEPTION("Unable to extend an incomplete data file at present.");
00217         }
00218 
00219         // Done!
00220         mIsInDefineMode = false;
00221         AdvanceAlongUnlimitedDimension();
00222     }
00223 }
00224 
00225 Hdf5DataWriter::~Hdf5DataWriter()
00226 {
00227     Close();
00228     
00229     if (mSinglePermutation)
00230     {
00231         MatDestroy(mSinglePermutation);
00232     }
00233     if (mDoublePermutation)
00234     {
00235         MatDestroy(mDoublePermutation);
00236     }
00237     if (mSingleIncompleteOutputMatrix)
00238     {
00239         MatDestroy(mSingleIncompleteOutputMatrix);
00240     }
00241     if (mDoubleIncompleteOutputMatrix)
00242     {
00243         MatDestroy(mDoubleIncompleteOutputMatrix);
00244     }
00245 }
00246 
00247 void Hdf5DataWriter::DefineFixedDimension(long dimensionSize)
00248 {
00249     if (!mIsInDefineMode)
00250     {
00251         EXCEPTION("Cannot define variables when not in Define mode");
00252     }
00253     if (dimensionSize < 1)
00254     {
00255         EXCEPTION("Fixed dimension must be at least 1 long");
00256     }
00257     if (mIsFixedDimensionSet)
00258     {
00259         EXCEPTION("Fixed dimension already set");
00260     }
00261 
00262     // Work out the ownership details
00263     mLo = mrVectorFactory.GetLow();
00264     mHi = mrVectorFactory.GetHigh();
00265     mNumberOwned = mrVectorFactory.GetLocalOwnership();
00266     mOffset = mLo;
00267     mFileFixedDimensionSize = dimensionSize;
00268     mDataFixedDimensionSize = dimensionSize;
00269     mIsFixedDimensionSet = true;
00270 }
00271 
00272 void Hdf5DataWriter::DefineFixedDimension(const std::vector<unsigned>& rNodesToOuput, long vecSize)
00273 {
00274     unsigned vector_size = rNodesToOuput.size();
00275 
00276     for (unsigned index=0; index < vector_size-1; index++)
00277     {
00278         if (rNodesToOuput[index] >= rNodesToOuput[index+1])
00279         {
00280             EXCEPTION("Input should be monotonic increasing");
00281         }
00282     }
00283 
00284     if ((int) rNodesToOuput.back() >= vecSize)
00285     {
00286         EXCEPTION("Vector size doesn't match nodes to output");
00287     }
00288 
00289     DefineFixedDimension(vecSize);
00290 
00291     mFileFixedDimensionSize = vector_size;
00292     mIsDataComplete = false;
00293     mIncompleteNodeIndices = rNodesToOuput;
00294     ComputeIncompleteOffset();
00295 }
00296 
00297 void Hdf5DataWriter::DefineFixedDimensionUsingMatrix(const std::vector<unsigned>& rNodesToOuput, long vecSize)
00298 {
00299     unsigned vector_size = rNodesToOuput.size();
00300 
00301     for (unsigned index=0; index < vector_size-1; index++)
00302     {
00303         if (rNodesToOuput[index] >= rNodesToOuput[index+1])
00304         {
00305             EXCEPTION("Input should be monotonic increasing");
00306         }
00307     }
00308 
00309     if ((int) rNodesToOuput.back() >= vecSize)
00310     {
00311         EXCEPTION("Vector size doesn't match nodes to output");
00312     }
00313     
00314     DefineFixedDimension(vecSize);
00315 
00316     mFileFixedDimensionSize = vector_size;
00317     mIsDataComplete = false;
00318     mIncompleteNodeIndices = rNodesToOuput;
00319     mUseMatrixForIncompleteData = true;
00320     ComputeIncompleteOffset();    
00321     
00322     //Make sure we've not done it already
00323     assert(mSingleIncompleteOutputMatrix == NULL);
00324     assert(mDoubleIncompleteOutputMatrix == NULL);
00325     PetscTools::SetupMat(mSingleIncompleteOutputMatrix,   mFileFixedDimensionSize,   mDataFixedDimensionSize, 2,  mNumberOwned,  mHi - mLo);
00326     PetscTools::SetupMat(mDoubleIncompleteOutputMatrix, 2*mFileFixedDimensionSize, 2*mDataFixedDimensionSize, 4,  2*mNumberOwned, 2*(mHi - mLo));
00327 
00328     //Only do local rows
00329     for (unsigned row_index = mOffset; row_index < mOffset + mNumberOwned; row_index++)
00330     {
00331         //Put zero on the diagonal
00332         MatSetValue(mSingleIncompleteOutputMatrix, row_index, row_index, 0.0, INSERT_VALUES);
00333         //Put one at (i,j)
00334         MatSetValue(mSingleIncompleteOutputMatrix, row_index, rNodesToOuput[row_index], 1.0, INSERT_VALUES);
00335         
00336         unsigned bi_index=2*row_index;
00337         unsigned perm_index=2*rNodesToOuput[row_index];
00338         //Put zeroes on the diagonal
00339         MatSetValue(mDoubleIncompleteOutputMatrix, bi_index, bi_index, 0.0, INSERT_VALUES);
00340         MatSetValue(mDoubleIncompleteOutputMatrix, bi_index+1, bi_index+1, 0.0, INSERT_VALUES);
00341         //Put ones at (i,j)
00342         MatSetValue(mDoubleIncompleteOutputMatrix, bi_index, perm_index, 1.0, INSERT_VALUES);
00343         MatSetValue(mDoubleIncompleteOutputMatrix, bi_index+1, perm_index+1, 1.0, INSERT_VALUES);
00344     }
00345     MatAssemblyBegin(mSingleIncompleteOutputMatrix, MAT_FINAL_ASSEMBLY);
00346     MatAssemblyBegin(mDoubleIncompleteOutputMatrix, MAT_FINAL_ASSEMBLY);
00347     MatAssemblyEnd(mSingleIncompleteOutputMatrix, MAT_FINAL_ASSEMBLY);
00348     MatAssemblyEnd(mDoubleIncompleteOutputMatrix, MAT_FINAL_ASSEMBLY);
00349     
00350 //    MatView(mSingleIncompleteOutputMatrix, PETSC_VIEWER_STDOUT_WORLD);
00351 
00352 }
00353 
00354 void Hdf5DataWriter::ComputeIncompleteOffset()
00355 {
00356     mOffset = 0;
00357     mNumberOwned = 0;
00358     for (unsigned i=0; i<mIncompleteNodeIndices.size(); i++)
00359     {
00360         if (mIncompleteNodeIndices[i] < mLo)
00361         {
00362             mOffset++;
00363         }
00364         else if (mIncompleteNodeIndices[i] < mHi)
00365         {
00366             mNumberOwned++;
00367         }
00368      }
00369 }
00370 
00371 int Hdf5DataWriter::DefineVariable(const std::string& rVariableName,
00372                                    const std::string& rVariableUnits)
00373 {
00374     if (!mIsInDefineMode)
00375     {
00376         EXCEPTION("Cannot define variables when not in Define mode");
00377     }
00378 
00379     CheckVariableName(rVariableName);
00380     CheckUnitsName(rVariableUnits);
00381 
00382     // Check for the variable being already defined
00383     for (unsigned index=0; index<mVariables.size(); index++)
00384     {
00385         if (mVariables[index].mVariableName == rVariableName)
00386         {
00387             EXCEPTION("Variable name already exists");
00388         }
00389     }
00390 
00391     DataWriterVariable new_variable;
00392     new_variable.mVariableName = rVariableName;
00393     new_variable.mVariableUnits = rVariableUnits;
00394     int variable_id;
00395 
00396     // Add the variable to the variable vector
00397     mVariables.push_back(new_variable);
00398 
00399     // Use the index of the variable vector as the variable ID.
00400     // This is ok since there is no way to remove variables.
00401     variable_id = mVariables.size() - 1;
00402 
00403     return variable_id;
00404 }
00405 
00406 void Hdf5DataWriter::CheckVariableName(const std::string& rName)
00407 {
00408     if (rName.length() == 0)
00409     {
00410         EXCEPTION("Variable name not allowed: may not be blank.");
00411     }
00412     CheckUnitsName(rName);
00413 }
00414 
00415 void Hdf5DataWriter::CheckUnitsName(const std::string& rName)
00416 {
00417     for (unsigned i=0; i<rName.length(); i++)
00418     {
00419         if (!isalnum(rName[i]) && !(rName[i]=='_'))
00420         {
00421             std::string error = "Variable name/units '" + rName + "' not allowed: may only contain alphanumeric characters or '_'.";
00422             EXCEPTION(error);
00423         }
00424     }
00425 }
00426 
00427 void Hdf5DataWriter::EndDefineMode()
00428 {
00429     //Check that at least one variable has been defined
00430     if (mVariables.size() < 1)
00431     {
00432         EXCEPTION("Cannot end define mode. No variables have been defined.");
00433     }
00434 
00435     //Check that a fixed dimension has been defined
00436     if (mIsFixedDimensionSet == false)
00437     {
00438         EXCEPTION("Cannot end define mode. One fixed dimension should be defined.");
00439     }
00440 
00441     OutputFileHandler output_file_handler(mDirectory, mCleanDirectory);
00442     std::string results_dir = output_file_handler.GetOutputDirectoryFullPath();
00443     std::string file_name = results_dir + mBaseName + ".h5";
00444 
00445     // Set up a property list saying how we'll open the file
00446     hid_t property_list_id = H5Pcreate(H5P_FILE_ACCESS);
00447     H5Pset_fapl_mpio(property_list_id, PETSC_COMM_WORLD, MPI_INFO_NULL);
00448 
00449     // Create a file (collectively) and free the property list
00450     mFileId = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, property_list_id);
00451     H5Pclose(property_list_id);
00452     if (mFileId < 0)
00453     {
00454         EXCEPTION("Hdf5DataWriter could not create " + file_name);
00455     }
00456     mIsInDefineMode = false;
00457 
00458     /*
00459      *  Create "Data" dataset
00460      */
00461     // Create the dataspace for the dataset.
00462     mDatasetDims[0] = 1; // While developing we got a non-documented "only the first dimension can be extendible" error.
00463     mDatasetDims[1] = mFileFixedDimensionSize;
00464     mDatasetDims[2] = mVariables.size();
00465 
00466     hsize_t* max_dims = NULL;
00467     hsize_t dataset_max_dims[DATASET_DIMS]; // dataset max dimensions
00468 
00469     hid_t cparms = H5P_DEFAULT;
00470 
00471     if (mIsUnlimitedDimensionSet)
00472     {
00473         dataset_max_dims[0] = H5S_UNLIMITED;
00474         dataset_max_dims[1] = mDatasetDims[1];
00475         dataset_max_dims[2] = mDatasetDims[2];
00476         max_dims = dataset_max_dims;
00477 
00478         hsize_t chunk_size;
00479         if (mUseOptimalChunkSizeAlgorithm)
00480         {
00481             // Modify dataset creation properties to enable chunking.
00482             // We don't want more than 100 chunks, as performance degrades significantly if there are too many,
00483             // where "too many" appears to be about 1000.
00484             // HDF5's caching won't apply if the chunks are too large, but this seems to have less of an impact.
00485             chunk_size = mEstimatedUnlimitedLength/100;
00486             if (chunk_size < 100)
00487             {
00488                 chunk_size = 100;
00489             }
00490         }
00491         else
00492         {
00493             chunk_size = mFixedChunkSize;
00494         }
00495         //If the size of a chunk in bytes is bigger than 4GB then there may be problems.
00496         // HDF5 1.6.x does not check for this - but there may be snags further down the line
00497         // HDF5 1.8.x does more error checking and produces std::cout errors and a file with no data in it
00498         if ( chunk_size * mDatasetDims[1] * mDatasetDims[2] >  (uint64_t)0xffffffff)
00499         {
00500             // Note: this exception can be avoided by altering the lines above
00501             // where chunk_size is set (at a loss of efficiency).
00502             mIsInDefineMode = true; //To stop things that would be created below from being deleted on Close()
00503             H5Fclose(mFileId); //This is the one thing which we have made
00504             EXCEPTION("HDF5 may be writing more than 4GB to disk at any time and would fail. It may be possible to tune the Chaste code to get around this");
00505         }
00506 
00507         hsize_t chunk_dims[DATASET_DIMS] = {chunk_size, mDatasetDims[1], mDatasetDims[2]};
00508         cparms = H5Pcreate (H5P_DATASET_CREATE);
00509         
00510         H5Pset_chunk( cparms, DATASET_DIMS, chunk_dims);
00511     }
00512 
00513     hid_t filespace = H5Screate_simple(DATASET_DIMS, mDatasetDims, max_dims);
00514 
00515     // Create the dataset and close filespace
00516     mDatasetId = H5Dcreate(mFileId, "Data", H5T_NATIVE_DOUBLE, filespace, cparms);
00517     H5Sclose(filespace);
00518 
00519     // Create dataspace for the name, unit attribute
00520     const unsigned MAX_STRING_SIZE = 100;
00521     hsize_t columns[1] = {mVariables.size()};
00522     hid_t colspace = H5Screate_simple(1, columns, NULL);
00523 
00524     // Create attribute for variable names
00525     char* col_data = (char*) malloc(mVariables.size() * sizeof(char) * MAX_STRING_SIZE);
00526 
00527     char* col_data_offset = col_data;
00528     for (unsigned var=0; var<mVariables.size(); var++)
00529     {
00530         std::string full_name = mVariables[var].mVariableName + "(" + mVariables[var].mVariableUnits + ")";
00531         strcpy (col_data_offset, full_name.c_str());
00532         col_data_offset += sizeof(char) * MAX_STRING_SIZE;
00533     }
00534 
00535     // Create the type 'string'
00536     hid_t string_type = H5Tcopy(H5T_C_S1);
00537     H5Tset_size(string_type, MAX_STRING_SIZE );
00538     hid_t attr = H5Acreate(mDatasetId, "Variable Details", string_type, colspace, H5P_DEFAULT);
00539 
00540     // Write to the attribute
00541     H5Awrite(attr, string_type, col_data);
00542 
00543     // Close dataspace & attribute
00544     free(col_data);
00545     H5Sclose(colspace);
00546     H5Aclose(attr);
00547 
00548     // Create "boolean" attribute telling the data to be incomplete or not
00549     columns[0] = 1;
00550     colspace = H5Screate_simple(1, columns, NULL);
00551     attr = H5Acreate(mDatasetId, "IsDataComplete", H5T_NATIVE_UINT, colspace, H5P_DEFAULT);
00552 
00553     // Write to the attribute - note that the native boolean is not predictable
00554     unsigned is_data_complete = mIsDataComplete ? 1 : 0;
00555     H5Awrite(attr, H5T_NATIVE_UINT, &is_data_complete);
00556 
00557     H5Sclose(colspace);
00558     H5Aclose(attr);
00559 
00560     if (!mIsDataComplete)
00561     {
00562         // We need to write a map
00563         // Create "unsigned" attribute with the map
00564         columns[0] = mFileFixedDimensionSize;
00565         colspace = H5Screate_simple(1, columns, NULL);
00566         attr = H5Acreate(mDatasetId, "NodeMap", H5T_NATIVE_UINT, colspace, H5P_DEFAULT);
00567 
00568         // Write to the attribute
00569         H5Awrite(attr, H5T_NATIVE_UINT, &mIncompleteNodeIndices[0]);
00570 
00571         H5Sclose(colspace);
00572         H5Aclose(attr);
00573     }
00574 
00575     /*
00576      *  Create "Time" dataset
00577      */
00578     if (mIsUnlimitedDimensionSet)
00579     {
00580         hsize_t time_dataset_dims[1] = {1};
00581         hsize_t time_dataset_max_dims[1] = {H5S_UNLIMITED};
00582 
00583         // Modify dataset creation properties to enable chunking.
00584         hsize_t time_chunk_dims[1] ={1};
00585         hid_t time_cparms = H5Pcreate (H5P_DATASET_CREATE);
00586         H5Pset_chunk( time_cparms, 1, time_chunk_dims);
00587 
00588         hid_t time_filespace = H5Screate_simple(1, time_dataset_dims, time_dataset_max_dims);
00589 
00590         // Create the dataset
00591         mTimeDatasetId = H5Dcreate(mFileId, mUnlimitedDimensionName.c_str(), H5T_NATIVE_DOUBLE, time_filespace, time_cparms);
00592 
00593         // Create the dataspace for the attribute
00594         hsize_t one = 1;
00595         hid_t one_column_space = H5Screate_simple(1, &one, NULL);
00596 
00597         // Create an attribute for the time unit
00598         hid_t unit_attr = H5Acreate(mTimeDatasetId, "Unit", string_type, one_column_space, H5P_DEFAULT);
00599 
00600         // Copy the unit to a string MAX_STRING_SIZE long and write it
00601         char unit_string[MAX_STRING_SIZE];
00602         strcpy(unit_string, mUnlimitedDimensionUnit.c_str());
00603         H5Awrite(unit_attr, string_type, unit_string);
00604 
00605         // Close the filespace
00606         H5Sclose(one_column_space);
00607         H5Aclose(unit_attr);
00608         H5Sclose(time_filespace);
00609     }
00610 
00611     /*
00612      * Create the provenance attribute
00613      */
00614     // Create a longer type for 'string'
00615     const unsigned MAX_PROVENANCE_STRING_SIZE = 1023;
00616     hid_t long_string_type = H5Tcopy(H5T_C_S1);
00617     H5Tset_size(long_string_type, MAX_PROVENANCE_STRING_SIZE );
00618     hsize_t prov_columns[1] = {1};
00619     hid_t provenance_space = H5Screate_simple(1, prov_columns, NULL);
00620     char* provenance_data = (char*) malloc(sizeof(char) * MAX_PROVENANCE_STRING_SIZE);
00621     assert( ChasteBuildInfo::GetProvenanceString().length() < MAX_PROVENANCE_STRING_SIZE);
00622     
00623     strcpy(provenance_data,  ChasteBuildInfo::GetProvenanceString().c_str());
00624     hid_t prov_attr = H5Acreate(mDatasetId, "Chaste Provenance", long_string_type, provenance_space, H5P_DEFAULT);
00625 
00626     // Write to the attribute
00627     H5Awrite(prov_attr, long_string_type, provenance_data);
00628 
00629     // Close dataspace & attribute
00630     free(provenance_data);
00631     H5Sclose(provenance_space);
00632     H5Aclose(prov_attr);
00633 
00634 }
00635 
00636 void Hdf5DataWriter::PutVector(int variableID, Vec petscVector)
00637 {
00638     if (mIsInDefineMode)
00639     {
00640         EXCEPTION("Cannot write data while in define mode.");
00641     }
00642 
00643     int vector_size;
00644     VecGetSize(petscVector, &vector_size);
00645 
00646     if ((unsigned) vector_size != mDataFixedDimensionSize)
00647     {
00648         EXCEPTION("Vector size doesn't match fixed dimension");
00649     }
00650 
00651     // Make sure that everything is actually extended to the correct dimension.
00652     PossiblyExtend();
00653     
00654     Vec output_petsc_vector;
00655     //Decide what to write 
00656     if (mSinglePermutation == NULL)
00657     {
00658         //No permutation - just write
00659         output_petsc_vector = petscVector;
00660     }
00661     else
00662     {
00663         assert(mIsDataComplete);
00664         //Make a vector with the same pattern (doesn't copy the data)
00665         VecDuplicate(petscVector, &output_petsc_vector);
00666         //Apply the permutation matrix
00667         MatMult(mSinglePermutation, petscVector, output_petsc_vector);
00668     }
00669     // Define a dataset in memory for this process
00670     hid_t memspace=0;
00671     if (mNumberOwned !=0)
00672     {
00673         hsize_t v_size[1] = {mNumberOwned};
00674         memspace = H5Screate_simple(1, v_size, NULL);
00675     }
00676     // Select hyperslab in the file
00677     hsize_t count[DATASET_DIMS] = {1, mNumberOwned, 1};
00678     hsize_t offset_dims[DATASET_DIMS] = {mCurrentTimeStep, mOffset, variableID};
00679     hid_t file_dataspace = H5Dget_space(mDatasetId);
00680 
00681     // Create property list for collective dataset
00682     hid_t property_list_id = H5Pcreate(H5P_DATASET_XFER);
00683     H5Pset_dxpl_mpio(property_list_id, H5FD_MPIO_COLLECTIVE);
00684 
00685     H5Sselect_hyperslab(file_dataspace, H5S_SELECT_SET, offset_dims, NULL, count, NULL);
00686 
00687     double* p_petsc_vector;
00688     VecGetArray(output_petsc_vector, &p_petsc_vector);
00689 
00690     if (mIsDataComplete)
00691     {
00692         H5Dwrite(mDatasetId, H5T_NATIVE_DOUBLE, memspace, file_dataspace, property_list_id, p_petsc_vector);
00693     }
00694     else
00695     {
00696         if (mUseMatrixForIncompleteData)
00697         {
00698             //Make a vector of the required size
00699             output_petsc_vector = PetscTools::CreateVec(mFileFixedDimensionSize, mNumberOwned);
00700           
00701             //Fill the vector by multiplying complete data by incomplete output matrix
00702             MatMult(mSingleIncompleteOutputMatrix, petscVector, output_petsc_vector);
00703             
00704             double* p_petsc_vector_incomplete;
00705             VecGetArray(output_petsc_vector, &p_petsc_vector_incomplete);
00706             H5Dwrite(mDatasetId, H5T_NATIVE_DOUBLE, memspace, file_dataspace, property_list_id, p_petsc_vector_incomplete);
00707         }
00708         else
00709         {
00710             // Make a local copy of the data you own
00711             double local_data[mNumberOwned];
00712             for (unsigned i=0; i<mNumberOwned; i++)
00713             {
00714                 local_data[i] = p_petsc_vector[ mIncompleteNodeIndices[mOffset+i]-mLo ];
00715     
00716             }
00717             H5Dwrite(mDatasetId, H5T_NATIVE_DOUBLE, memspace, file_dataspace, property_list_id, local_data);
00718         }
00719     }
00720 
00721     VecRestoreArray(output_petsc_vector, &p_petsc_vector);
00722 
00723     H5Pclose(property_list_id);
00724     H5Sclose(file_dataspace);
00725     if (mNumberOwned !=0)
00726     {
00727         H5Sclose(memspace);
00728     }
00729     
00730     if (petscVector != output_petsc_vector)
00731     {
00732         //Free local vector
00733         VecDestroy(output_petsc_vector);
00734     }
00735 }
00736 
00737 void Hdf5DataWriter::PutStripedVector(std::vector<int> variableIDs, Vec petscVector)
00738 {
00739     if (mIsInDefineMode)
00740     {
00741         EXCEPTION("Cannot write data while in define mode.");
00742     }
00743 
00744     if (variableIDs.size() <= 1)
00745     {
00746         EXCEPTION("The PutStripedVector method requires at least two variables ID. If only one is needed, use PutVector method instead");
00747     }
00748 
00749     const int NUM_STRIPES=variableIDs.size();
00750 
00751     int firstVariableID=variableIDs[0];
00752 
00753     // Currently the method only works with consecutive columns, can be extended if needed
00754     for (unsigned i = 1; i < variableIDs.size(); i++)
00755     {
00756         if (variableIDs[i]-variableIDs[i-1] != 1)
00757         {
00758             EXCEPTION("Columns should be consecutive. Try reordering them.");
00759         }
00760     }
00761 
00762     int vector_size;
00763     VecGetSize(petscVector, &vector_size);
00764 
00765     if ((unsigned) vector_size != NUM_STRIPES*mDataFixedDimensionSize)
00766     {
00767         EXCEPTION("Vector size doesn't match fixed dimension");
00768     }
00769 
00770     // Make sure that everything is actually extended to the correct dimension
00771     PossiblyExtend();
00772 
00773     Vec output_petsc_vector;
00774     //Decide what to write 
00775     if (mDoublePermutation == NULL)
00776     {
00777         //No permutation - just write
00778         output_petsc_vector = petscVector;
00779     }
00780     else
00781     {
00782         assert(mIsDataComplete);
00783         //Make a vector with the same pattern (doesn't copy the data)
00784         VecDuplicate(petscVector, &output_petsc_vector);
00785         //Apply the permutation matrix
00786         MatMult(mDoublePermutation, petscVector, output_petsc_vector);
00787     }
00788     // Define a dataset in memory for this process
00789     hid_t memspace=0;
00790     if (mNumberOwned !=0)
00791     {
00792         hsize_t v_size[1] = {mNumberOwned*NUM_STRIPES};
00793         memspace = H5Screate_simple(1, v_size, NULL);
00794     }
00795 
00796     // Select hyperslab in the file
00797     hsize_t start[DATASET_DIMS] = {mCurrentTimeStep, mOffset, firstVariableID};
00798     hsize_t stride[DATASET_DIMS] = {1, 1, 1};//we are imposing contiguous variables, hence the stride is 1 (3rd component)
00799     hsize_t block_size[DATASET_DIMS] = {1, mNumberOwned, 1};
00800     hsize_t number_blocks[DATASET_DIMS] = {1, 1, NUM_STRIPES};
00801 
00802     hid_t hyperslab_space = H5Dget_space(mDatasetId);
00803     H5Sselect_hyperslab(hyperslab_space, H5S_SELECT_SET, start, stride, number_blocks, block_size);
00804 
00805     // Create property list for collective dataset write, and write! Finally.
00806     hid_t property_list_id = H5Pcreate(H5P_DATASET_XFER);
00807     H5Pset_dxpl_mpio(property_list_id, H5FD_MPIO_COLLECTIVE);
00808 
00809     double* p_petsc_vector;
00810     VecGetArray(output_petsc_vector, &p_petsc_vector);
00811 
00812     if (mIsDataComplete)
00813     {
00814         H5Dwrite(mDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, p_petsc_vector);
00815     }
00816     else
00817     {
00818         if(variableIDs.size() < 3) //incomplete data and striped vector is supported only for NUM_STRIPES=2...for the moment
00819         {
00820             if (mUseMatrixForIncompleteData)
00821             {
00822                 //Make a vector of the required size
00823                 output_petsc_vector = PetscTools::CreateVec(2*mFileFixedDimensionSize, 2*mNumberOwned);
00824 
00825                 //Fill the vector by multiplying complete data by incomplete output matrix
00826                 MatMult(mDoubleIncompleteOutputMatrix, petscVector, output_petsc_vector);
00827 
00828                 double* p_petsc_vector_incomplete;
00829                 VecGetArray(output_petsc_vector, &p_petsc_vector_incomplete);
00830 
00831                 H5Dwrite(mDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, p_petsc_vector_incomplete);
00832             }
00833             else
00834             {
00835                 // Make a local copy of the data you own
00836                 double local_data[mNumberOwned*NUM_STRIPES];
00837                 for (unsigned i=0; i<mNumberOwned; i++)
00838                 {
00839                     unsigned local_node_number = mIncompleteNodeIndices[mOffset+i] - mLo;
00840                     local_data[NUM_STRIPES*i]   = p_petsc_vector[ local_node_number*NUM_STRIPES ];
00841                     local_data[NUM_STRIPES*i+1] = p_petsc_vector[ local_node_number*NUM_STRIPES + 1];
00842                 }
00843 
00844                 H5Dwrite(mDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, local_data);
00845             }
00846         }
00847         else
00848         {
00849             EXCEPTION("The PutStripedVector functionality for incomplete data is supported for only 2 stripes");
00850         }
00851     }
00852 
00853     VecRestoreArray(output_petsc_vector, &p_petsc_vector);
00854 
00855     H5Sclose(hyperslab_space);
00856     if (mNumberOwned != 0)
00857     {
00858         H5Sclose(memspace);
00859     }
00860     H5Pclose(property_list_id);
00861     
00862     if (petscVector != output_petsc_vector)
00863     {
00864         //Free local vector
00865         VecDestroy(output_petsc_vector);
00866     }
00867 }
00868 
00869 void Hdf5DataWriter::PutUnlimitedVariable(double value)
00870 {
00871     if (mIsInDefineMode)
00872     {
00873         EXCEPTION("Cannot write data while in define mode.");
00874     }
00875 
00876     if (!mIsUnlimitedDimensionSet)
00877     {
00878         EXCEPTION("PutUnlimitedVariable() called but no unlimited dimension has been set");
00879     }
00880 
00881     // Make sure that everything is actually extended to the correct dimension.
00882     PossiblyExtend();
00883 
00884     // This data is only written by the master
00885     if (!PetscTools::AmMaster())
00886     {
00887         return;
00888     }
00889 
00890     hsize_t size[1] = {1};
00891     hid_t memspace = H5Screate_simple(1, size, NULL);
00892 
00893     // Select hyperslab in the file.
00894     hsize_t count[1] = {1};
00895     hsize_t offset[1] = {mCurrentTimeStep};
00896     hid_t hyperslab_space = H5Dget_space(mTimeDatasetId);
00897     H5Sselect_hyperslab(hyperslab_space, H5S_SELECT_SET, offset, NULL, count, NULL);
00898 
00899     H5Dwrite(mTimeDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, H5P_DEFAULT, &value);
00900 
00901     H5Sclose(hyperslab_space);
00902     H5Sclose(memspace);
00903 }
00904 
00905 void Hdf5DataWriter::Close()
00906 {
00907     if (mIsInDefineMode)
00908     {
00909         return; // Nothing to do...
00910     }
00911     
00912     H5Dclose(mDatasetId);
00913     if (mIsUnlimitedDimensionSet)
00914     {
00915         H5Dclose(mTimeDatasetId);
00916     }
00917     H5Fclose(mFileId);
00918     
00919     // Cope with being called twice (e.g. if a user calls Close then the destructor)
00920     mIsInDefineMode = true;
00921 }
00922 
00923 void Hdf5DataWriter::DefineUnlimitedDimension(const std::string& rVariableName,
00924                                               const std::string& rVariableUnits,
00925                                               unsigned estimatedLength)
00926 {
00927     if (mIsUnlimitedDimensionSet)
00928     {
00929         EXCEPTION("Unlimited dimension already set. Cannot be defined twice");
00930     }
00931 
00932     if (!mIsInDefineMode)
00933     {
00934         EXCEPTION("Cannot define variables when not in Define mode");
00935     }
00936 
00937     mIsUnlimitedDimensionSet = true;
00938     mUnlimitedDimensionName = rVariableName;
00939     mUnlimitedDimensionUnit = rVariableUnits;
00940     mEstimatedUnlimitedLength = estimatedLength;
00941 }
00942 
00943 void Hdf5DataWriter::AdvanceAlongUnlimitedDimension()
00944 {
00945     if (!mIsUnlimitedDimensionSet)
00946     {
00947         EXCEPTION("Trying to advance along an unlimited dimension without having defined any");
00948     }
00949 
00950     mCurrentTimeStep++;
00951 
00952     /*
00953      *  Extend the dataset
00954      */
00955     // If the user provided an estimate for the length of the
00956     // unlimited dimension, allocate all that space.
00957     if ( mEstimatedUnlimitedLength > mDatasetDims[0] )
00958     {
00959         mDatasetDims[0] = mEstimatedUnlimitedLength;
00960         mNeedExtend = true;
00961     }
00962 
00963     // If you are beyond the user estimate increment step by step
00964     if ( mCurrentTimeStep >= (long) mEstimatedUnlimitedLength )
00965     {
00966         mDatasetDims[0]++;
00967         mNeedExtend = true;
00968     }
00969 
00970 
00971 }
00972 
00973 void Hdf5DataWriter::PossiblyExtend()
00974 {
00975     if (mNeedExtend)
00976     {
00977         H5Dextend (mDatasetId, mDatasetDims);
00978         H5Dextend (mTimeDatasetId, mDatasetDims);
00979     }
00980     mNeedExtend = false;
00981 }
00982 
00983 int Hdf5DataWriter::GetVariableByName(const std::string& rVariableName)
00984 {
00985     int id = -1;
00986     // Check for the variable name in the existing variables
00987     for (unsigned index=0; index<mVariables.size(); index++)
00988     {
00989         if (mVariables[index].mVariableName == rVariableName)
00990         {
00991             id = index;
00992             break;
00993         }
00994     }
00995     if (id == -1)
00996     {
00997         EXCEPTION("Variable does not exist in hdf5 definitions.");
00998     }
00999     return id;
01000 }
01001 
01002 bool Hdf5DataWriter::ApplyPermutation(const std::vector<unsigned>& rPermutation)
01003 {
01004     if (!mIsInDefineMode)
01005     {
01006         EXCEPTION("Cannot define permutation when not in Define mode");
01007     }
01008     
01009     if (rPermutation.empty())
01010     {
01011         return false;
01012     }
01013     
01014     if (rPermutation.size() != mFileFixedDimensionSize || 
01015         rPermutation.size() != mDataFixedDimensionSize)
01016     {
01017         EXCEPTION("Permutation doesn't match the expected problem size");
01018     }
01019     //Permutation checker
01020     std::set<unsigned> permutation_pigeon_hole;
01021     //Fill up the pigeon holes
01022     bool identity_map=true;
01023     for (unsigned i=0; i<mDataFixedDimensionSize; i++)
01024     {
01025         permutation_pigeon_hole.insert(rPermutation[i]);
01026         if (identity_map && i != rPermutation[i])
01027         {
01028             identity_map=false;
01029         }
01030     }
01031     if (identity_map)
01032     {
01033         //Do nothing
01034         return false;
01035     }
01036     /* Pigeon-hole principal says that each index appears exactly once
01037      * so if any don't appear then either one appears twice or something out of
01038      * scope has appeared. 
01039      */
01040     for (unsigned i=0; i<mDataFixedDimensionSize; i++)
01041     {
01042         if (permutation_pigeon_hole.count(i) != 1u)
01043         {
01044             EXCEPTION("Permutation vector doesn't contain a valid permutation");  
01045         }
01046     }
01047     //Make sure we've not done it already
01048     assert(mSinglePermutation == NULL);
01049     assert(mDoublePermutation == NULL);
01050     PetscTools::SetupMat(mSinglePermutation,   mDataFixedDimensionSize,   mDataFixedDimensionSize, 2, mHi - mLo, mHi - mLo);
01051     PetscTools::SetupMat(mDoublePermutation, 2*mDataFixedDimensionSize, 2*mDataFixedDimensionSize, 4, 2*(mHi - mLo), 2*(mHi - mLo));
01052 #if (PETSC_VERSION_MAJOR == 3) //PETSc 3.x.x
01053     MatSetOption(mSinglePermutation, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE); 
01054     MatSetOption(mDoublePermutation, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE); 
01055 #else
01056     MatSetOption(mSinglePermutation, MAT_IGNORE_OFF_PROC_ENTRIES); 
01057     MatSetOption(mDoublePermutation, MAT_IGNORE_OFF_PROC_ENTRIES); 
01058 #endif
01059     //Only do local rows
01060     for (unsigned row_index=mLo; row_index<mHi; row_index++)
01061     {
01062         //Put zero on the diagonal
01063         MatSetValue(mSinglePermutation, row_index, row_index, 0.0, INSERT_VALUES);
01064         //Put one at (i,j)
01065         MatSetValue(mSinglePermutation, row_index, rPermutation[row_index], 1.0, INSERT_VALUES);
01066         
01067         unsigned bi_index=2*row_index;
01068         unsigned perm_index=2*rPermutation[row_index];
01069         //Put zeroes on the diagonal
01070         MatSetValue(mDoublePermutation, bi_index, bi_index, 0.0, INSERT_VALUES);
01071         MatSetValue(mDoublePermutation, bi_index+1, bi_index+1, 0.0, INSERT_VALUES);
01072         //Put ones at (i,j)
01073         MatSetValue(mDoublePermutation, bi_index, perm_index, 1.0, INSERT_VALUES);
01074         MatSetValue(mDoublePermutation, bi_index+1, perm_index+1, 1.0, INSERT_VALUES);
01075     }
01076     MatAssemblyBegin(mSinglePermutation, MAT_FINAL_ASSEMBLY);
01077     MatAssemblyBegin(mDoublePermutation, MAT_FINAL_ASSEMBLY);
01078     MatAssemblyEnd(mSinglePermutation, MAT_FINAL_ASSEMBLY);
01079     MatAssemblyEnd(mDoublePermutation, MAT_FINAL_ASSEMBLY);  
01080     return true;  
01081 }
01082 
01083 void Hdf5DataWriter::SetFixedChunkSize(unsigned chunkSize)
01084 {
01085     assert(mIsInDefineMode);
01086 
01087     mUseOptimalChunkSizeAlgorithm = false;
01088     mFixedChunkSize = chunkSize;
01089 }
01090 

Generated on Tue May 31 14:31:42 2011 for Chaste by  doxygen 1.5.5