Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/feature/h5-read-reservoir-pressu…
Browse files Browse the repository at this point in the history 
…re' into feature/Small-fixes-to-file-handling-files
  • Loading branch information
bellout committed Jan 4, 2017
2 parents e99373c + 178458b commit fb1095c
Show file tree
Hide file tree
Showing 3 changed files with 199 additions and 28 deletions.
116 changes: 105 additions & 11 deletions FieldOpt/Hdf5SummaryReader/hdf5_summary_reader.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,27 +8,89 @@ Hdf5SummaryReader::Hdf5SummaryReader(const std::string file_path)
DATASET_NAME_TIMES("TIMES"),
GROUP_NAME_FLOW_TRANSPORT("FLOW_TRANSPORT"),
DATASET_NAME_ACTIVE_CELLS("ACTIVE_CELLS"),
DATASET_NAME_WELL_STATES("WELL_STATES")
DATASET_NAME_WELL_STATES("WELL_STATES"),
DATASET_NAME_PRESSURE("PTZ")
{
readTimeVector(file_path);
readActiveCells(file_path);
readWellStates(file_path);
readActiveCells(file_path);
readReservoirPressure(file_path);
}

void Hdf5SummaryReader::readActiveCells(std::string file_path) {
void Hdf5SummaryReader::readReservoirPressure(std::string file_path) {

// Read the file
H5File file(file_path, H5F_ACC_RDONLY);
Group group = Group(file.openGroup(GROUP_NAME_FLOW_TRANSPORT));
DataSet dataset = DataSet(group.openDataSet(DATASET_NAME_ACTIVE_CELLS));
DataSet dataset = DataSet(group.openDataSet(DATASET_NAME_PRESSURE));

DataSpace dataspace = dataset.getSpace();
hsize_t dims[2];
hsize_t dims[3];

auto rank = dataspace.getSimpleExtentDims(dims, NULL);
// Uncomment to debug:
// std::cout << "dataset rank = " << rank << ", dimensions "
// << (unsigned long)(dims[0]) << " x "
// << (unsigned long)(dims[1]) << std::endl;
// << (unsigned long)(dims[1]) << " x "
// << (unsigned long)(dims[2]) << std::endl;

// Define hyperslab
hsize_t count[3] = {dims[0], 1, dims[2]};
hsize_t offset[3] = {0, 0, 0};
dataspace.selectHyperslab(H5S_SELECT_SET, count, offset);

// Size of selected column + define memory space
hsize_t col_sz[2] = {dims[0], dims[2]};
DataSpace mspace(2, col_sz);

std::vector<double> vector;
vector.resize(dims[0] * dims[2]);
dataset.read(vector.data(), PredType::NATIVE_DOUBLE, mspace, dataspace);

// Reorder pressure vector
// Note:
// Data/time component ordering inside vector:
// Example: pressure vector with 5 cells over 8 time steps:
// size of vector = ncells (x ndata_cols=1) x ntime_steps
//
// cell 1 cell 2 cell 3 cell 4 cell 5
//time |--------|--------|--------|--------|--------|
//steps: 12345678 12345678 12345678 12345678 12345678

int ncells = (int)dims[0];
int ndata_cols = (int)dims[1];
int ntime_steps = (int)dims[2];

for (int tt = 0; tt < ntime_steps; ++tt) {
std::vector<double> sub;
for (int cc = 0; cc < ncells; ++cc) {
sub.push_back(vector[ cc * ntime_steps + tt ]);
}
pressure_.push_back(sub);

// Uncomment to debug:
// std::cout << "pressure_[tt=" << tt << "].size() = " << pressure_[tt].size() << std::endl;
}

// Uncomment to debug:
// for (int rr = 0; rr < 10; ++rr) { // pressure_.size()
// std::cout << "pressure = " << pressure_[0][rr] << std::endl;
// }

}

void Hdf5SummaryReader::readActiveCells(std::string file_path) {

// Read the file
H5File file(file_path, H5F_ACC_RDONLY);
Group group = Group(file.openGroup(GROUP_NAME_FLOW_TRANSPORT));
DataSet dataset = DataSet(group.openDataSet(DATASET_NAME_ACTIVE_CELLS));

DataSpace dataspace = dataset.getSpace();
hsize_t dims[2];

// rank variable can be used for debug
auto rank = dataspace.getSimpleExtentDims(dims, NULL);

// Define hyperslab
hsize_t count[2] = { dims[0], 1 };
Expand All @@ -39,10 +101,13 @@ void Hdf5SummaryReader::readActiveCells(std::string file_path) {
hsize_t col_sz[1] = { dims[0] };
DataSpace mspace( 1, col_sz );

// Read and reorder vector
std::vector<int> vector;
vector.resize(dims[0]);
dataset.read(vector.data(), PredType::NATIVE_INT, mspace, dataspace);
active_cells_ = vector;

cells_all_vector_ = vector;
cells_find_statuses(cells_all_vector_);
}

void Hdf5SummaryReader::readTimeVector(std::string file_path) {
Expand All @@ -60,10 +125,22 @@ void Hdf5SummaryReader::readTimeVector(std::string file_path) {
hsize_t dims[2];
dataspace.getSimpleExtentDims(dims, NULL);

// Check size of time vector is > 1, otherwise tricky errors
// occurs further downstream in readWellStates()
if ( (int)dims[0] < 2) {
throw std::runtime_error("TIME vector has only one component! "
"Check your simulation H5 output.");
}

std::vector<double> vector;
vector.resize(dims[0]);
dataset.read(vector.data(), PredType::NATIVE_DOUBLE);
times_ = vector;

// Uncomment to debug:
// std::cout << "size of time vector = "
// << vector.size() << std::endl;

}

void Hdf5SummaryReader::readWellStates(std::string file_path) {
Expand All @@ -76,8 +153,8 @@ void Hdf5SummaryReader::readWellStates(std::string file_path) {
hsize_t dims[2];
dataspace.getSimpleExtentDims(dims, NULL);

nwells_ = dims[0];
ntimes_ = dims[1];
nwells_ = (int)dims[0];
ntimes_ = (int)dims[1];

CompType ctype(sizeof(wstype_t));
auto double_type = PredType::NATIVE_DOUBLE;
Expand Down Expand Up @@ -127,9 +204,9 @@ void Hdf5SummaryReader::parseWsVector(std::vector<wstype_t> &wsvec) {
Hdf5SummaryReader::well_data Hdf5SummaryReader::parseWellState(std::vector<wstype_t> &ws, int wnr) {
int first = wnr*ntimes_;
int last = first + ntimes_;
int nperfs = ws[first].vAverageDensity.size();
int nperfs = (int)ws[first].vAverageDensity.size();
auto state = Hdf5SummaryReader::well_data(ntimes_, nperfs);
state.nphases = ws[first].vPhaseRates.size() / nperfs;
state.nphases = (int)ws[first].vPhaseRates.size() / nperfs;
int t = 0;
for (int i = first; i < last; ++i) { // Well data at each time step
state.well_types[t] = ws[i].vIntData[0];
Expand Down Expand Up @@ -259,6 +336,23 @@ int Hdf5SummaryReader::number_of_phases() const {
return nphases_;
}

int Hdf5SummaryReader::cells_find_statuses(std::vector<int> &cells_all_vector_) {

for (int i = 0; i < cells_all_vector_.size(); ++i) {
if (cells_all_vector_[i] < 0){
cells_inactive_.push_back(cells_all_vector_[i]);
cells_inactive_idx_.push_back(i);
}else{
cells_active_.push_back(cells_all_vector_[i]);
cells_active_idx_.push_back(i);
}
}

cells_total_num_ = (int)cells_all_vector_.size();
cells_num_active_ = (int)cells_active_.size();
cells_num_inactive_ = (int)cells_inactive_.size();
}

std::vector<double> Hdf5SummaryReader::field_cumulative_oil_production_sc() const {
std::vector<double> sum(ntimes_, 0.0);
for (int w = 0; w < nwells_; ++w) {
Expand Down
71 changes: 58 additions & 13 deletions FieldOpt/Hdf5SummaryReader/hdf5_summary_reader.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,8 @@ class Hdf5SummaryReader {
/*!
* Read the HDF5 summary file written by AD-GPRS at the specified path.
*
* \note The provided path is not checked by this method, and should therefore be checked before
* invoking this.
* \note The provided path is not checked by this method, and should therefore be
* checked before invoking this.
* @param file_path Path to a .H5 summary file.
* @return A Hdf5SummaryReader object containing information from the summary.
*/
Expand All @@ -39,9 +39,34 @@ class Hdf5SummaryReader {
const std::vector<double> &times_steps() const { return times_; }

/*!
* Get the vector containing defining active cells.
* Get reservoir pressure vector.
*/
const std::vector<int> &active_cells() const { return active_cells_; }
std::vector<std::vector<double>> reservoir_pressure() const { return pressure_; }

/*!
* Return vector of active grid cells.
*/
const std::vector<int> &cells_active() { return cells_active_; }

/*!
* Return vector of active grid cell indices.
*/
const std::vector<int> &cells_active_idx() { return cells_active_idx_; }

/*!
* Get total number of grid cells in model.
*/
int cells_total_num() const { return cells_total_num_; }

/*!
* Get total number of active grid cells in model.
*/
int cells_num_active() const { return cells_num_active_; }

/*!
* Get total number of inactive grid cells in model.
*/
int cells_num_inactive() const { return cells_num_inactive_; }

/*!
* Get the number of wells found in the summary.
Expand Down Expand Up @@ -185,16 +210,19 @@ class Hdf5SummaryReader {
std::vector<double> field_cumulative_gas_injection_sc() const;



private:
const H5std_string GROUP_NAME_RESTART; //!< The name of the restart group in the HDF5 file.
const H5std_string DATASET_NAME_TIMES; //!< The name of the dataset containing the time step vector in the HDF5 file.
const H5std_string GROUP_NAME_FLOW_TRANSPORT; //!< The name of the flow transport group in the HDF5 file.
const H5std_string DATASET_NAME_WELL_STATES; //!< The name of the dataset containing the well states in the HDF5 file.
const H5std_string DATASET_NAME_ACTIVE_CELLS; //!< The name of the dataset containing active cells vector.
const H5std_string DATASET_NAME_PRESSURE; //!< The name of the dataset containing pressure and component data.

/*!
* The wstype_t datatype represents the datatype in which well states are stored in the HDF5 file. Each element
* in the dataset contains information about a well and its perforations at a specific time step.
* The wstype_t datatype represents the datatype in which well states are stored in the HDF5 file.
* Each element in the dataset contains information about a well and its perforations at a specific
* time step.
*/
typedef struct wstype_t {
std::vector<double> vPressures;
Expand All @@ -214,8 +242,8 @@ class Hdf5SummaryReader {
} wstype_t;

/*!
* The perforation_data struct holds vectors containing rate, pressure, temperature and density values at
* all time steps for a specific perforation in a well.
* The perforation_data struct holds vectors containing rate, pressure, temperature and density
* values at all time steps for a specific perforation in a well.
*/
struct perforation_data {
perforation_data(){}
Expand All @@ -229,8 +257,9 @@ class Hdf5SummaryReader {
};

/*!
* The well_data struct holds information about a specific well, as well as a vector of perforation_data for
* all of its perforations, and vectors containing rate and pressure values at all time steps.
* The well_data struct holds information about a specific well, as well as a vector of
* perforation_data for all of its perforations, and vectors containing rate and pressure
* values at all time steps.
*/
struct well_data {
well_data(){}
Expand All @@ -249,16 +278,32 @@ class Hdf5SummaryReader {
};

void readTimeVector(std::string file_path); //!< Read the time vector from the HDF5 summary file.
void readActiveCells(std::string file_path); //!< Read vector defining which cells are active from the HDF5 summary file.
void readWellStates(std::string file_path); //!< Read all well state information from the HDF5 summary file.
void parseWsVector(std::vector<wstype_t> &wsvec); //!< Populate well_states_ by creating well_data and perforation_data objects from the wstype_t vector.
well_data parseWellState(std::vector<wstype_t> &ws, int wnr); //!< Parse the states for a single well and create a well_data object.

void readActiveCells(std::string file_path); //!< Read vector defining which cells are active from the HDF5 summary file.
void readReservoirPressure(std::string file_path); //!< Read reservoir cell pressures for each time step.

/*!
* variables containing information about the reservoir cell ensemble,
* e.g., number of active cells, corresponding active cell indices, etc.
*/
std::vector<int> cells_all_vector_; //!< Vector def. status for all cells (size equal to total number of cells).
std::vector<int> cells_active_, cells_inactive_;
std::vector<int> cells_active_idx_, cells_inactive_idx_;
int cells_total_num_; //!< Total number of grid cells in model.
int cells_num_active_; //!< Total number of active grid cells in model.
int cells_num_inactive_; //!< Total number of inactive grid cells in model.

int cells_find_statuses(std::vector<int> &cells_all_vector_); //!< Fills inactive/active cells numbers and indices

int nwells_; //!< Number of wells in summary.
int ntimes_; //!< Number of time steps in the summary.
int nphases_; //!< Number of phases in the model
std::vector<int> active_cells_; //!< Vector defining active cells (size equal to total number of cells) .
int nphases_; //!< Number of phases in the model.

std::vector<double> times_; //!< Vector containing all time steps.
std::vector<std::vector<double>> pressure_; //!< Vector containing reservoir pressures.

/*!
* Calculate the cumulative using the composite trapezoidal rule.
Expand Down
40 changes: 36 additions & 4 deletions FieldOpt/Hdf5SummaryReader/tests/test_hdf5_summary_reader.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include "../hdf5_summary_reader.h"
#include <gtest/gtest.h>
#include <hdf5_summary_reader.h>

namespace {
class Hdf5SummaryReaderTest : public ::testing::Test {
Expand Down Expand Up @@ -35,12 +36,43 @@ namespace {

TEST_F(Hdf5SummaryReaderTest, ActiveCellVector) {
auto reader = Hdf5SummaryReader(file_path);
std::vector<int> active_cells = reader.active_cells();

for (int i = 0; i < active_cells.size(); ++i) {
// std::cout << active_cells[i] << std::endl;
EXPECT_EQ(active_cells[i], i);
auto cells_total_num_ = reader.cells_total_num();
auto cells_num_active_ = reader.cells_num_active();
auto cells_num_inactive_ = reader.cells_num_inactive();

auto cells_active_ = reader.cells_active();
auto cells_active_idx_ = reader.cells_active_idx();

// Test based on all cells in test reservoir (5spot) being active
EXPECT_EQ(cells_total_num_, 3600);
EXPECT_EQ(cells_num_active_, cells_total_num_);
EXPECT_EQ(cells_num_inactive_, 0);

for (int i = 0; i < cells_num_active_; ++i) {
// std::cout << cells_num_active_[i] << std::endl;
EXPECT_EQ(cells_active_[i], i);
EXPECT_EQ(cells_active_idx_[i], i);
}
}

TEST_F(Hdf5SummaryReaderTest, readReservoirPressure) {
auto reader = Hdf5SummaryReader(file_path);
auto pressure = reader.reservoir_pressure();

std::vector<double> def_pressure = {
370.555603, 370.55451121327735, 366.80490627652864,
281.07546615360917, 213.85427339240036, 179.59397120055957,
154.40334284097705, 151.17045510746573
};

// Test against first component of each pressure vector
// to confirm we're reading correct column and that the
// read vector has been correctly resized
for (int i = 0; i < pressure.size(); ++i) {
EXPECT_EQ(def_pressure[i], pressure[i][0]);
}

}

TEST_F(Hdf5SummaryReaderTest, IntegerData) {
Expand Down

0 comments on commit fb1095c

Please sign in to comment.