diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp index 8aaab2711..b672cc30b 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp @@ -27,38 +27,17 @@ namespace Opm const UnstructuredGrid& grid, bool init_rock) { - if (init_rock){ - rock_.init(deck, grid); - } - pvt_.init(deck, 0); - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(deck, grid, 0); - - if (pvt_.numPhases() != satprops_->numPhases()) { - OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); - } + init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids, + grid.dimensions, init_rock); } + BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, const UnstructuredGrid& grid, bool init_rock) { - if (init_rock){ - rock_.init(newParserDeck, grid); - } - pvt_.init(newParserDeck, /*numSamples=*/200); - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(newParserDeck, grid, /*numSamples=*/200); - - if (pvt_.numPhases() != satprops_->numPhases()) { - OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); - } + init(newParserDeck, grid.number_of_cells, grid.global_cell, grid.cartdims, + grid.cell_centroids, grid.dimensions, init_rock); } BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, @@ -66,63 +45,8 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock) { - if(init_rock){ - rock_.init(deck, grid); - } - - const int pvt_samples = param.getDefault("pvt_tab_size", 0); - pvt_.init(deck, pvt_samples); - - // Unfortunate lack of pointer smartness here... - const int sat_samples = param.getDefault("sat_tab_size", 0); - std::string threephase_model = param.getDefault("threephase_model", "simple"); - if (deck.hasField("ENDSCALE") && threephase_model != "gwseg") { - OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'gwseg' model only."); - } - if (sat_samples > 1) { - if (threephase_model == "stone2") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(deck, grid, sat_samples); - } else if (threephase_model == "simple") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(deck, grid, sat_samples); - } else if (threephase_model == "gwseg") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(deck, grid, sat_samples); - } else { - OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model); - } - } else { - if (threephase_model == "stone2") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(deck, grid, sat_samples); - } else if (threephase_model == "simple") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(deck, grid, sat_samples); - } else if (threephase_model == "gwseg") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(deck, grid, sat_samples); - } else { - OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model); - } - } - - if (pvt_.numPhases() != satprops_->numPhases()) { - OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); - } + init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids, + grid.dimensions, param, init_rock); } BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, @@ -130,63 +54,8 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock) { - if(init_rock){ - rock_.init(newParserDeck, grid); - } - - const int pvt_samples = param.getDefault("pvt_tab_size", 200); - pvt_.init(newParserDeck, pvt_samples); - - // Unfortunate lack of pointer smartness here... - const int sat_samples = param.getDefault("sat_tab_size", 200); - std::string threephase_model = param.getDefault("threephase_model", "simple"); - if (newParserDeck->hasKeyword("ENDSCALE") && threephase_model != "simple") { - OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'simple' model only."); - } - if (sat_samples > 1) { - if (threephase_model == "stone2") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(newParserDeck, grid, sat_samples); - } else if (threephase_model == "simple") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(newParserDeck, grid, sat_samples); - } else if (threephase_model == "gwseg") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(newParserDeck, grid, sat_samples); - } else { - OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model); - } - } else { - if (threephase_model == "stone2") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(newParserDeck, grid, sat_samples); - } else if (threephase_model == "simple") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(newParserDeck, grid, sat_samples); - } else if (threephase_model == "gwseg") { - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(newParserDeck, grid, sat_samples); - } else { - OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model); - } - } - - if (pvt_.numPhases() != satprops_->numPhases()) { - OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); - } + init(newParserDeck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids, + grid.dimensions, param, init_rock); } BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck() diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index 5c74a76a6..e521e28aa 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -90,6 +90,45 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock=true); + template + BlackoilPropertiesFromDeck(const EclipseGridParser& deck, + int number_of_cells, + const int* global_cell, + const int* cart_dims, + T begin_cell_centroids, + int dimension, + bool init_rock=true); + + + template + BlackoilPropertiesFromDeck(const EclipseGridParser& deck, + int number_of_cells, + const int* global_cell, + const int* cart_dims, + T begin_cell_centroids, + int dimension, + const parameter::ParameterGroup& param, + bool init_rock=true); + + template + BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, + int number_of_cells, + const int* global_cell, + const int* cart_dims, + T begin_cell_centroids, + int dimension, + bool init_rock=true); + + template + BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, + int number_of_cells, + const int* global_cell, + const int* cart_dims, + T begin_cell_centroids, + int dimension, + const parameter::ParameterGroup& param, + bool init_rock=true); + /// Destructor. virtual ~BlackoilPropertiesFromDeck(); @@ -211,6 +250,41 @@ namespace Opm double* smax) const; private: + template + void init(const EclipseGridParser& deck, + int number_of_cells, + const int* global_cell, + const int* cart_dims, + T begin_cell_centroids, + int dimension, + bool init_rock); + template + void init(const EclipseGridParser& deck, + int number_of_cells, + const int* global_cell, + const int* cart_dims, + T begin_cell_centroids, + int dimension, + const parameter::ParameterGroup& param, + bool init_rock); + + template + void init(Opm::DeckConstPtr newParserDeck, + int number_of_cells, + const int* global_cell, + const int* cart_dims, + T begin_cell_centroids, + int dimension, + bool init_rock); + template + void init(Opm::DeckConstPtr newParserDeck, + int number_of_cells, + const int* global_cell, + const int* cart_dims, + T begin_cell_centroids, + int dimension, + const parameter::ParameterGroup& param, + bool init_rock); RockFromDeck rock_; BlackoilPvtProperties pvt_; std::unique_ptr satprops_; @@ -224,5 +298,6 @@ namespace Opm } // namespace Opm +#include "BlackoilPropertiesFromDeck_impl.hpp" #endif // OPM_BLACKOILPROPERTIESFROMDECK_HEADER_INCLUDED diff --git a/opm/core/props/rock/RockFromDeck.cpp b/opm/core/props/rock/RockFromDeck.cpp index 1058fa997..e34d784da 100644 --- a/opm/core/props/rock/RockFromDeck.cpp +++ b/opm/core/props/rock/RockFromDeck.cpp @@ -65,45 +65,58 @@ namespace Opm void RockFromDeck::init(const EclipseGridParser& deck, const UnstructuredGrid& grid) { - assignPorosity(deck, grid); - permfield_valid_.assign(grid.number_of_cells, false); + init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims); + } + + /// Initialize from deck and cell mapping. + /// \param deck Deck input parser + /// \param number_of_cells The number of cells in the grid. + /// \param global_cell The mapping fom local to global cell indices. + /// global_cell[i] is the corresponding global index of i. + /// \param cart_dims The size of the underlying cartesian grid. + void RockFromDeck::init(const EclipseGridParser& deck, + int number_of_cells, const int* global_cell, + const int* cart_dims) + { + assignPorosity(deck, number_of_cells, global_cell); + permfield_valid_.assign(number_of_cells, false); const double perm_threshold = 0.0; // Maybe turn into parameter? - assignPermeability(deck, grid, perm_threshold); + assignPermeability(deck, number_of_cells, global_cell, cart_dims, perm_threshold); } void RockFromDeck::init(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid) + int number_of_cells, const int* global_cell, + const int* cart_dims) { - assignPorosity(newParserDeck, grid); - permfield_valid_.assign(grid.number_of_cells, false); + assignPorosity(newParserDeck, number_of_cells, global_cell); + permfield_valid_.assign(number_of_cells, false); const double perm_threshold = 0.0; // Maybe turn into parameter? - assignPermeability(newParserDeck, grid, perm_threshold); + assignPermeability(newParserDeck, number_of_cells, global_cell, cart_dims, + perm_threshold); } void RockFromDeck::assignPorosity(const EclipseGridParser& parser, - const UnstructuredGrid& grid) + int number_of_cells, const int* global_cell) { - porosity_.assign(grid.number_of_cells, 1.0); - const int* gc = grid.global_cell; + porosity_.assign(number_of_cells, 1.0); if (parser.hasField("PORO")) { const std::vector& poro = parser.getFloatingPointValue("PORO"); for (int c = 0; c < int(porosity_.size()); ++c) { - const int deck_pos = (gc == NULL) ? c : gc[c]; + const int deck_pos = (global_cell == NULL) ? c : global_cell[c]; porosity_[c] = poro[deck_pos]; } } } void RockFromDeck::assignPorosity(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid) + int number_of_cells, const int* global_cell) { - porosity_.assign(grid.number_of_cells, 1.0); - const int* gc = grid.global_cell; + porosity_.assign(number_of_cells, 1.0); if (newParserDeck->hasKeyword("PORO")) { const std::vector& poro = newParserDeck->getKeyword("PORO")->getSIDoubleData(); for (int c = 0; c < int(porosity_.size()); ++c) { - const int deck_pos = (gc == NULL) ? c : gc[c]; + const int deck_pos = (global_cell == NULL) ? c : global_cell[c]; assert(0 <= c && c < (int) porosity_.size()); assert(0 <= deck_pos && deck_pos < (int) poro.size()); porosity_[c] = poro[deck_pos]; @@ -113,16 +126,17 @@ namespace Opm void RockFromDeck::assignPermeability(const EclipseGridParser& parser, - const UnstructuredGrid& grid, + int number_of_cells, + const int* global_cell, + const int* cartdims, double perm_threshold) { const int dim = 3; - const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2]; - const int nc = grid.number_of_cells; + const int num_global_cells = cartdims[0]*cartdims[1]*cartdims[2]; assert(num_global_cells > 0); - permeability_.assign(dim * dim * nc, 0.0); + permeability_.assign(dim * dim * number_of_cells, 0.0); std::vector*> tensor; tensor.reserve(10); @@ -144,13 +158,12 @@ namespace Opm // chosen) default value... // if (tensor.size() > 1) { - const int* gc = grid.global_cell; int off = 0; - for (int c = 0; c < nc; ++c, off += dim*dim) { + for (int c = 0; c < number_of_cells; ++c, off += dim*dim) { // SharedPermTensor K(dim, dim, &permeability_[off]); int kix = 0; - const int glob = (gc == NULL) ? c : gc[c]; + const int glob = (global_cell == NULL) ? c : global_cell[c]; for (int i = 0; i < dim; ++i) { for (int j = 0; j < dim; ++j, ++kix) { @@ -167,12 +180,14 @@ namespace Opm } void RockFromDeck::assignPermeability(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid, + int number_of_cells, + const int* global_cell, + const int* cartdims, double perm_threshold) { const int dim = 3; - const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2]; - const int nc = grid.number_of_cells; + const int num_global_cells = cartdims[0]*cartdims[1]*cartdims[2]; + const int nc = number_of_cells; assert(num_global_cells > 0); @@ -198,7 +213,7 @@ namespace Opm // chosen) default value... // if (tensor.size() > 1) { - const int* gc = grid.global_cell; + const int* gc = global_cell; int off = 0; for (int c = 0; c < nc; ++c, off += dim*dim) { diff --git a/opm/core/props/rock/RockFromDeck.hpp b/opm/core/props/rock/RockFromDeck.hpp index b5e0fb4e0..5d8c44aef 100644 --- a/opm/core/props/rock/RockFromDeck.hpp +++ b/opm/core/props/rock/RockFromDeck.hpp @@ -46,13 +46,24 @@ namespace Opm void init(const EclipseGridParser& deck, const UnstructuredGrid& grid); - /// Initialize from deck and grid. - /// \param newParserDeck Deck produced by the opm-parser code - /// \param grid Grid to which property object applies, needed for the - /// mapping from cell indices (typically from a processed grid) - /// to logical cartesian indices consistent with the deck. + /// Initialize from deck and cell mapping. + /// \param deck Deck input parser + /// \param number_of_cells The number of cells in the grid. + /// \param global_cell The mapping fom local to global cell indices. + /// global_cell[i] is the corresponding global index of i. + /// \param cart_dims The size of the underlying cartesian grid. + void init(const EclipseGridParser& deck, + int number_of_cells, const int* global_cell, + const int* cart_dims); + /// Initialize from deck and cell mapping. + /// \param newParserDeck Deck produced by the opm-parser code + /// \param number_of_cells The number of cells in the grid. + /// \param global_cell The mapping fom local to global cell indices. + /// global_cell[i] is the corresponding global index of i. + /// \param cart_dims The size of the underlying cartesian grid. void init(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid); + int number_of_cells, const int* global_cell, + const int* cart_dims); /// \return D, the number of spatial dimensions. Always 3 for deck input. int numDimensions() const @@ -82,14 +93,20 @@ namespace Opm private: void assignPorosity(const EclipseGridParser& parser, - const UnstructuredGrid& grid); + int number_of_cells, + const int* global_cell); void assignPorosity(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid); + int number_of_cells, + const int* global_cell); void assignPermeability(const EclipseGridParser& parser, - const UnstructuredGrid& grid, + int number_of_cells, + const int* global_cell, + const int* cart_dims, const double perm_threshold); void assignPermeability(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid, + int number_of_cells, + const int* global_cell, + const int* cart_dims, double perm_threshold); std::vector porosity_; diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp index 757b1481a..bf52b3061 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp @@ -64,7 +64,29 @@ namespace Opm const int samples); /// Initialize from deck and grid. - /// \param[in] deck Deck input parser + /// \param[in] deck Deck input parser + /// \param[in] number_of_cells The number of cells of the grid to which property + /// object applies, needed for the + /// mapping from cell indices (typically from a processed + /// grid) to logical cartesian indices consistent with the + /// deck. + /// \param[in] global_cell The mapping from local cell indices of the grid to + /// global cell indices used in the deck. + /// \param[in] begin_cell_centroids Pointer to the first cell_centroid of the grid. + /// \param[in] dimensions The dimensions of the grid. + /// \param[in] samples Number of uniform sample points for saturation tables. + /// \tparam T The iterator Type for the cell centroids. + /// NOTE: samples will only be used with the SatFuncSetUniform template argument. + template + void init(const EclipseGridParser& deck, + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroids, + int dimensions, + const int samples); + + /// Initialize from deck and grid. + /// \param[in] newParserDeck Deck input parser /// \param[in] grid Grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. @@ -74,6 +96,27 @@ namespace Opm const UnstructuredGrid& grid, const int samples); + /// Initialize from deck and grid. + /// \param[in] newParserDeck Deck input parser + /// \param[in] number_of_cells The number of cells of the grid to which property + /// object applies, needed for the + /// mapping from cell indices (typically from a processed + /// grid) to logical cartesian indices consistent with the + /// deck. + /// \param[in] global_cell The mapping from local cell indices of the grid to + /// global cell indices used in the deck. + /// \param[in] begin_cell_centroids Pointer to the first cell_centroid of the grid. + /// \param[in] dimensions The dimensions of the grid. + /// \param[in] samples Number of uniform sample points for saturation tables. + /// NOTE: samples will only be used with the SatFuncSetUniform template argument. + template + void init(Opm::DeckConstPtr newParserDeck, + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroids, + int dimensions, + const int samples); + /// \return P, the number of phases. int numPhases() const; @@ -138,20 +181,45 @@ namespace Opm typedef SatFuncSet Funcs; const Funcs& funcForCell(const int cell) const; + + template void initEPS(const EclipseGridParser& deck, - const UnstructuredGrid& grid); + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroids, + int dimensions); + template void initEPSHyst(const EclipseGridParser& deck, - const UnstructuredGrid& grid); + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroids, + int dimensions); + template void initEPSKey(const EclipseGridParser& deck, - const UnstructuredGrid& grid, + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroids, + int dimensions, const std::string& keyword, std::vector& scaleparam); + template void initEPS(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid); + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroids, + int dimensions); + template void initEPSHyst(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid); + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroids, + int dimensions); + template void initEPSKey(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid, + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroids, + int dimensions, const std::string& keyword, std::vector& scaleparam); void initEPSParam(const int cell, diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp index c85cef78c..54d716da6 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -51,6 +52,20 @@ namespace Opm void SaturationPropsFromDeck::init(const EclipseGridParser& deck, const UnstructuredGrid& grid, const int samples) + { + init(deck, grid.number_of_cells, grid.global_cell, grid.cell_centroids, + grid.dimensions, samples); + } + + /// Initialize from deck. + template + template< class T> + void SaturationPropsFromDeck::init(const EclipseGridParser& deck, + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroids, + int dimensions, + const int samples) { phase_usage_ = phaseUsageFromDeck(deck); @@ -66,11 +81,9 @@ namespace Opm if (deck.hasField("SATNUM")) { const std::vector& satnum = deck.getIntegerValue("SATNUM"); satfuncs_expected = *std::max_element(satnum.begin(), satnum.end()); - const int num_cells = grid.number_of_cells; - cell_to_func_.resize(num_cells); - const int* gc = grid.global_cell; - for (int cell = 0; cell < num_cells; ++cell) { - const int deck_pos = (gc == NULL) ? cell : gc[cell]; + cell_to_func_.resize(number_of_cells); + for (int cell = 0; cell < number_of_cells; ++cell) { + const int deck_pos = (global_cell == NULL) ? cell : global_cell[cell]; cell_to_func_[cell] = satnum[deck_pos] - 1; } } @@ -124,8 +137,7 @@ namespace Opm } } do_eps_ = true; - - initEPS(deck, grid); + initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions); // For now, a primitive detection of hysteresis. TODO: SATOPTS HYSTER/ and EHYSTR do_hyst_ = deck.hasField("ISWL") || deck.hasField("ISWU") || deck.hasField("ISWCR") || deck.hasField("ISGL") || @@ -137,7 +149,8 @@ namespace Opm deck.hasField("IKRGR") || deck.hasField("IKRORW") || deck.hasField("IKRORG") ) { OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Currently hysteresis and relperm value scaling can not be combined."); } - initEPSHyst(deck, grid); + initEPSHyst(deck, number_of_cells, global_cell, begin_cell_centroids, + dimensions); } //OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Under construction ..."); @@ -150,6 +163,21 @@ namespace Opm void SaturationPropsFromDeck::init(Opm::DeckConstPtr newParserDeck, const UnstructuredGrid& grid, const int samples) + { + this->template init(newParserDeck, grid.number_of_cells, + grid.global_cell, grid.cell_centroids, + grid.dimensions, samples); + } + + /// Initialize from deck. + template + template + void SaturationPropsFromDeck::init(Opm::DeckConstPtr newParserDeck, + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroids, + int dimensions, + const int samples) { phase_usage_ = phaseUsageFromDeck(newParserDeck); @@ -165,9 +193,9 @@ namespace Opm if (newParserDeck->hasKeyword("SATNUM")) { const std::vector& satnum = newParserDeck->getKeyword("SATNUM")->getIntData(); satfuncs_expected = *std::max_element(satnum.begin(), satnum.end()); - const int num_cells = grid.number_of_cells; + const int num_cells = number_of_cells; cell_to_func_.resize(num_cells); - const int* gc = grid.global_cell; + const int* gc = global_cell; for (int cell = 0; cell < num_cells; ++cell) { const int newParserDeck_pos = (gc == NULL) ? cell : gc[cell]; cell_to_func_[cell] = satnum[newParserDeck_pos] - 1; @@ -225,7 +253,8 @@ namespace Opm } do_eps_ = true; - initEPS(newParserDeck, grid); + initEPS(newParserDeck, number_of_cells, global_cell, begin_cell_centroids, + dimensions); // For now, a primitive detection of hysteresis. TODO: SATOPTS HYSTER/ and EHYSTR do_hyst_ = @@ -257,7 +286,8 @@ namespace Opm "Currently hysteresis and relperm value scaling " "cannot be combined."); } - initEPSHyst(newParserDeck, grid); + initEPSHyst(newParserDeck, number_of_cells, global_cell, begin_cell_centroids, + dimensions); } } } @@ -441,31 +471,51 @@ namespace Opm return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]]; } - // Initialize saturation scaling parameters + + // Initialize saturation scaling parameter template + template void SaturationPropsFromDeck::initEPS(const EclipseGridParser& deck, - const UnstructuredGrid& grid) + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroid, + int dimensions) { std::vector swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr; std::vector krw, krg, kro, krwr, krgr, krorw, krorg; // Initialize saturation scaling parameter - initEPSKey(deck, grid, std::string("SWL"), swl); - initEPSKey(deck, grid, std::string("SWU"), swu); - initEPSKey(deck, grid, std::string("SWCR"), swcr); - initEPSKey(deck, grid, std::string("SGL"), sgl); - initEPSKey(deck, grid, std::string("SGU"), sgu); - initEPSKey(deck, grid, std::string("SGCR"), sgcr); - initEPSKey(deck, grid, std::string("SOWCR"), sowcr); - initEPSKey(deck, grid, std::string("SOGCR"), sogcr); - initEPSKey(deck, grid, std::string("KRW"), krw); - initEPSKey(deck, grid, std::string("KRG"), krg); - initEPSKey(deck, grid, std::string("KRO"), kro); - initEPSKey(deck, grid, std::string("KRWR"), krwr); - initEPSKey(deck, grid, std::string("KRGR"), krgr); - initEPSKey(deck, grid, std::string("KRORW"), krorw); - initEPSKey(deck, grid, std::string("KRORG"), krorg); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SWL"), swl); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SWU"), swu); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SWCR"), swcr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SGL"), sgl); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SGU"), sgu); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SGCR"), sgcr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SOWCR"), sowcr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SOGCR"), sogcr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRW"), krw); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRG"), krg); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRO"), kro); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRWR"), krwr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRGR"), krgr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRORW"), krorw); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRORG"), krorg); - eps_transf_.resize(grid.number_of_cells); + eps_transf_.resize(number_of_cells); const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; @@ -473,7 +523,7 @@ namespace Opm const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; - for (int cell = 0; cell < grid.number_of_cells; ++cell) { + for (int cell = 0; cell < number_of_cells; ++cell) { if (oilWater) { // ### krw initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], @@ -507,29 +557,48 @@ namespace Opm // Initialize saturation scaling parameters template + template void SaturationPropsFromDeck::initEPS(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid) + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroid, + int dimensions) { std::vector swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr; std::vector krw, krg, kro, krwr, krgr, krorw, krorg; // Initialize saturation scaling parameter - initEPSKey(newParserDeck, grid, std::string("SWL"), swl); - initEPSKey(newParserDeck, grid, std::string("SWU"), swu); - initEPSKey(newParserDeck, grid, std::string("SWCR"), swcr); - initEPSKey(newParserDeck, grid, std::string("SGL"), sgl); - initEPSKey(newParserDeck, grid, std::string("SGU"), sgu); - initEPSKey(newParserDeck, grid, std::string("SGCR"), sgcr); - initEPSKey(newParserDeck, grid, std::string("SOWCR"), sowcr); - initEPSKey(newParserDeck, grid, std::string("SOGCR"), sogcr); - initEPSKey(newParserDeck, grid, std::string("KRW"), krw); - initEPSKey(newParserDeck, grid, std::string("KRG"), krg); - initEPSKey(newParserDeck, grid, std::string("KRO"), kro); - initEPSKey(newParserDeck, grid, std::string("KRWR"), krwr); - initEPSKey(newParserDeck, grid, std::string("KRGR"), krgr); - initEPSKey(newParserDeck, grid, std::string("KRORW"), krorw); - initEPSKey(newParserDeck, grid, std::string("KRORG"), krorg); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SWL"), swl); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SWU"), swu); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SWCR"), swcr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SGL"), sgl); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SGU"), sgu); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SGCR"), sgcr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SOWCR"), sowcr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("SOGCR"), sogcr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRW"), krw); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRG"), krg); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRO"), kro); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRWR"), krwr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRGR"), krgr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRORW"), krorw); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("KRORG"), krorg); - eps_transf_.resize(grid.number_of_cells); + eps_transf_.resize(number_of_cells); const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; @@ -537,7 +606,7 @@ namespace Opm const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; - for (int cell = 0; cell < grid.number_of_cells; ++cell) { + for (int cell = 0; cell < number_of_cells; ++cell) { if (oilWater) { // ### krw initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], @@ -571,31 +640,50 @@ namespace Opm // Initialize hysteresis saturation scaling parameters template + template void SaturationPropsFromDeck::initEPSHyst(const EclipseGridParser& deck, - const UnstructuredGrid& grid) + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroid, + int dimensions) { std::vector iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr; std::vector ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg; // Initialize hysteresis saturation scaling parameters - initEPSKey(deck, grid, std::string("ISWL"), iswl); - initEPSKey(deck, grid, std::string("ISWU"), iswu); - initEPSKey(deck, grid, std::string("ISWCR"), iswcr); - initEPSKey(deck, grid, std::string("ISGL"), isgl); - initEPSKey(deck, grid, std::string("ISGU"), isgu); - initEPSKey(deck, grid, std::string("ISGCR"), isgcr); - initEPSKey(deck, grid, std::string("ISOWCR"), isowcr); - initEPSKey(deck, grid, std::string("ISOGCR"), isogcr); - initEPSKey(deck, grid, std::string("IKRW"), ikrw); - initEPSKey(deck, grid, std::string("IKRG"), ikrg); - initEPSKey(deck, grid, std::string("IKRO"), ikro); - initEPSKey(deck, grid, std::string("IKRWR"), ikrwr); - initEPSKey(deck, grid, std::string("IKRGR"), ikrgr); - initEPSKey(deck, grid, std::string("IKRORW"), ikrorw); - initEPSKey(deck, grid, std::string("IKRORG"), ikrorg); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISWL"), iswl); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISWU"), iswu); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISWCR"), iswcr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISGL"), isgl); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISGU"), isgu); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISGCR"), isgcr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISOWCR"), isowcr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISOGCR"), isogcr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRW"), ikrw); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRG"), ikrg); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRO"), ikro); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRWR"), ikrwr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRGR"), ikrgr); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRORW"), ikrorw); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRORG"), ikrorg); - eps_transf_hyst_.resize(grid.number_of_cells); - sat_hyst_.resize(grid.number_of_cells); + eps_transf_hyst_.resize(number_of_cells); + sat_hyst_.resize(number_of_cells); const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; @@ -603,7 +691,7 @@ namespace Opm const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; - for (int cell = 0; cell < grid.number_of_cells; ++cell) { + for (int cell = 0; cell < number_of_cells; ++cell) { if (oilWater) { // ### krw initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], @@ -637,30 +725,49 @@ namespace Opm // Initialize hysteresis saturation scaling parameters template + template void SaturationPropsFromDeck::initEPSHyst(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid) + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroid, + int dimensions) { std::vector iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr; std::vector ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg; // Initialize hysteresis saturation scaling parameters - initEPSKey(newParserDeck, grid, std::string("ISWL"), iswl); - initEPSKey(newParserDeck, grid, std::string("ISWU"), iswu); - initEPSKey(newParserDeck, grid, std::string("ISWCR"), iswcr); - initEPSKey(newParserDeck, grid, std::string("ISGL"), isgl); - initEPSKey(newParserDeck, grid, std::string("ISGU"), isgu); - initEPSKey(newParserDeck, grid, std::string("ISGCR"), isgcr); - initEPSKey(newParserDeck, grid, std::string("ISOWCR"), isowcr); - initEPSKey(newParserDeck, grid, std::string("ISOGCR"), isogcr); - initEPSKey(newParserDeck, grid, std::string("IKRW"), ikrw); - initEPSKey(newParserDeck, grid, std::string("IKRG"), ikrg); - initEPSKey(newParserDeck, grid, std::string("IKRO"), ikro); - initEPSKey(newParserDeck, grid, std::string("IKRWR"), ikrwr); - initEPSKey(newParserDeck, grid, std::string("IKRGR"), ikrgr); - initEPSKey(newParserDeck, grid, std::string("IKRORW"), ikrorw); - initEPSKey(newParserDeck, grid, std::string("IKRORG"), ikrorg); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISWL"), iswl); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISWU"), iswu); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISWCR"), iswcr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISGL"), isgl); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISGU"), isgu); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISGCR"), isgcr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISOWCR"), isowcr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("ISOGCR"), isogcr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRW"), ikrw); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRG"), ikrg); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRO"), ikro); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRWR"), ikrwr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRGR"), ikrgr); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRORW"), ikrorw); + initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IKRORG"), ikrorg); - eps_transf_hyst_.resize(grid.number_of_cells); - sat_hyst_.resize(grid.number_of_cells); + eps_transf_hyst_.resize(number_of_cells); + sat_hyst_.resize(number_of_cells); const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; @@ -668,7 +775,7 @@ namespace Opm const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; - for (int cell = 0; cell < grid.number_of_cells; ++cell) { + for (int cell = 0; cell < number_of_cells; ++cell) { if (oilWater) { // ### krw initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], @@ -702,8 +809,12 @@ namespace Opm // Initialize saturation scaling parameter template + template void SaturationPropsFromDeck::initEPSKey(const EclipseGridParser& deck, - const UnstructuredGrid& grid, + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroid, + int dimensions, const std::string& keyword, std::vector& scaleparam) { @@ -724,57 +835,57 @@ namespace Opm if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) { if (useAqua && (useKeyword || deck.getENPTVD().mask_[0])) { itab = 1; - scaleparam.resize(grid.number_of_cells); - for (int i=0; i& val = deck.getFloatingPointValue(keyword); for (int c = 0; c < int(scaleparam.size()); ++c) { - const int deck_pos = (gc == NULL) ? c : gc[c]; + const int deck_pos = (global_cell == NULL) ? c : global_cell[c]; scaleparam[c] = val[deck_pos]; } } else { @@ -858,14 +968,15 @@ namespace Opm deck.getENPTVD().write(std::cout); else deck.getENKRVD().write(std::cout); - const double* cc = grid.cell_centroids; - const int dim = grid.dimensions; - for (int cell = 0; cell < grid.number_of_cells; ++cell) { + for (int cell = 0; cell < number_of_cells; ++cell) { int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell]; if (table[itab][jtab][0] != -1.0) { std::vector& depth = table[0][jtab]; std::vector& val = table[itab][jtab]; - double zc = cc[dim*cell+dim-1]; + double zc = UgGridHelpers + ::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, + dimensions), + dimensions-1); if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval scaleparam[cell] = linearInterpolation(depth, val, zc); } @@ -876,8 +987,12 @@ namespace Opm // Initialize saturation scaling parameter template + template void SaturationPropsFromDeck::initEPSKey(Opm::DeckConstPtr newParserDeck, - const UnstructuredGrid& grid, + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroid, + int dimensions, const std::string& keyword, std::vector& scaleparam) { @@ -898,57 +1013,57 @@ namespace Opm if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) { if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 0))) { itab = 1; - scaleparam.resize(grid.number_of_cells); - for (int i=0; i& val = newParserDeck->getKeyword(keyword)->getSIDoubleData(); for (int c = 0; c < int(scaleparam.size()); ++c) { const int deck_pos = (gc == NULL) ? c : gc[c]; @@ -1065,9 +1180,9 @@ namespace Opm else newParserDeck.getENKRVD().write(std::cout); */ - const double* cc = grid.cell_centroids; - const int dim = grid.dimensions; - for (int cell = 0; cell < grid.number_of_cells; ++cell) { + const double* cc = begin_cell_centroid; + const int dim = dimensions; + for (int cell = 0; cell < number_of_cells; ++cell) { int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell]; if (table[itab][jtab][0] != -1.0) { std::vector& depth = table[0][jtab]; diff --git a/opm/core/simulator/BlackoilState.cpp b/opm/core/simulator/BlackoilState.cpp index 0b864a17d..64f921bd2 100644 --- a/opm/core/simulator/BlackoilState.cpp +++ b/opm/core/simulator/BlackoilState.cpp @@ -4,13 +4,18 @@ using namespace Opm; void -BlackoilState::init(const UnstructuredGrid& g, int num_phases) { - SimulatorState::init(g, num_phases); - gor_.resize(g.number_of_cells, 0.) ; - rv_.resize(g.number_of_cells,0.); +BlackoilState::init(int number_of_cells, int number_of_phases, int num_phases) { + SimulatorState::init(number_of_cells, number_of_phases, num_phases); + gor_.resize(number_of_cells, 0.) ; + rv_.resize(number_of_cells,0.); // surfvol_ intentionally empty, left to initBlackoilSurfvol } +void +BlackoilState::init(const UnstructuredGrid& g, int num_phases) +{ + init(g.number_of_cells, g.number_of_faces, num_phases); +} /// Set the first saturation to either its min or max value in /// the indicated cells. The second saturation value s2 is set /// to (1.0 - s1) for each cell. Any further saturation values diff --git a/opm/core/simulator/BlackoilState.hpp b/opm/core/simulator/BlackoilState.hpp index c6fb7779e..0d55f4004 100644 --- a/opm/core/simulator/BlackoilState.hpp +++ b/opm/core/simulator/BlackoilState.hpp @@ -32,7 +32,9 @@ namespace Opm class BlackoilState : public SimulatorState { public: - virtual void init(const UnstructuredGrid& g, int num_phases); + virtual void init(const UnstructuredGrid& grid, int num_phases); + + virtual void init(int number_of_cells, int number_of_faces, int num_phases); /// Set the first saturation to either its min or max value in /// the indicated cells. The second saturation value s2 is set diff --git a/opm/core/simulator/initState.hpp b/opm/core/simulator/initState.hpp index 0cbaf769b..6bdbbef5d 100644 --- a/opm/core/simulator/initState.hpp +++ b/opm/core/simulator/initState.hpp @@ -65,6 +65,44 @@ namespace Opm const double gravity, State& state); + /// Initialize a two-phase state from parameters. + /// The following parameters are accepted (defaults): + /// - convection_testcase (false) -- Water in the 'left' part of the grid. + /// - ref_pressure (100) -- Initial pressure in bar for all cells + /// (if convection_testcase is true), + /// or pressure at woc depth. + /// - segregation_testcase (false) -- Water above the woc instead of below. + /// - water_oil_contact (none) -- Depth of water-oil contact (woc). + /// - init_saturation (none) -- Initial water saturation for all cells. + /// + /// If convection_testcase is true, the saturation is initialised + /// as indicated, and pressure is initialised to a constant value + /// ('ref_pressure'). + /// If segregation_testcase is true, the saturation is initialised + /// as indicated, and pressure is initialised hydrostatically. + /// Otherwise we have 3 cases: + /// 1. If 'water_oil_contact' is given, saturation is initialised + /// accordingly. + /// 2. If 'water_oil_contact' is not given, but 'init_saturation' + /// is given, water saturation is set to that value everywhere. + /// 3. If neither are given, water saturation is set to minimum. + /// + /// In all three cases, pressure is initialised hydrostatically. + /// In case 2) and 3), the depth of the first cell is used as reference depth. + template + void initStateBasic(int number_of_cells, + const int* global_cell, + const int* cartdims, + int number_of_faces, + FaceCells face_cells, + FCI begin_face_centroids, + CCI begin_cell_centroids, + int dimensions, + const IncompPropertiesInterface& props, + const parameter::ParameterGroup& param, + const double gravity, + State& state); + /// Initialize a blackoil state from parameters. /// The following parameters are accepted (defaults): /// - convection_testcase (false) -- Water in the 'left' part of the grid. @@ -88,6 +126,36 @@ namespace Opm const double gravity, State& state); + /// Initialize a blackoil state from parameters. + /// The following parameters are accepted (defaults): + /// - convection_testcase (false) -- Water in the 'left' part of the grid. + /// - ref_pressure (100) -- Initial pressure in bar for all cells + /// (if convection_testcase is true), + /// or pressure at woc depth. + /// - water_oil_contact (none) -- Depth of water-oil contact (woc). + /// If convection_testcase is true, the saturation is initialised + /// as indicated, and pressure is initialised to a constant value + /// ('ref_pressure'). + /// Otherwise we have 2 cases: + /// 1. If 'water_oil_contact' is given, saturation is initialised + /// accordingly. + /// 2. Water saturation is set to minimum. + /// In both cases, pressure is initialised hydrostatically. + /// In case 2., the depth of the first cell is used as reference depth. + template + void initStateBasic(int number_of_cells, + const int* global_cell, + const int* cartdims, + int number_of_faces, + FaceCells face_cells, + FCI begin_face_centroids, + CCI begin_cell_centroids, + int dimensions, + const BlackoilPropertiesInterface& props, + const parameter::ParameterGroup& param, + const double gravity, + State& state); + /// Initialize a two-phase state from input deck. /// If EQUIL is present: /// - saturation is set according to the water-oil contact, @@ -102,6 +170,27 @@ namespace Opm const double gravity, State& state); + /// Initialize a two-phase state from input deck. + /// If EQUIL is present: + /// - saturation is set according to the water-oil contact, + /// - pressure is set to hydrostatic equilibrium. + /// Otherwise: + /// - saturation is set according to SWAT, + /// - pressure is set according to PRESSURE. + template + void initStateFromDeck(int number_of_cells, + const int* global_cell, + const int* cartdims, + int number_of_faces, + FaceCells face_cells, + FCI begin_face_centroids, + CCI begin_cell_centroids, + int dimensions, + const Props& props, + const EclipseGridParser& deck, + const double gravity, + State& state); + /// Initialize a two-phase water-oil blackoil state from input deck. /// If EQUIL is present: /// - saturation is set according to the water-oil contact, @@ -117,7 +206,26 @@ namespace Opm const double gravity, State& state); - + /// Initialize a two-phase water-oil blackoil state from input deck. + /// If EQUIL is present: + /// - saturation is set according to the water-oil contact, + /// - pressure is set to hydrostatic equilibrium. + /// Otherwise: + /// - saturation is set according to SWAT, + /// - pressure is set according to PRESSURE. + /// In addition, this function sets surfacevol. + template + void initBlackoilStateFromDeck(int number_of_cells, + const int* global_cell, + int number_of_faces, + FaceCells face_cells, + FCI begin_face_centroids, + CCI begin_cell_centroids, + int dimensions, + const Props& props, + const EclipseGridParser& deck, + const double gravity, + State& state); } // namespace Opm #include diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index 18372573e..32758f0e3 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -49,17 +50,21 @@ namespace Opm // Find the cells that are below and above a depth. // TODO: add 'anitialiasing', obtaining a more precise split // by f. ex. subdividing cells cut by the split depths. - void cellsBelowAbove(const UnstructuredGrid& grid, + template + void cellsBelowAbove(int number_of_cells, + T begin_cell_centroids, + int dimension, const double depth, std::vector& below, std::vector& above) { - const int num_cells = grid.number_of_cells; - below.reserve(num_cells); - above.reserve(num_cells); - const int dim = grid.dimensions; - for (int c = 0; c < num_cells; ++c) { - const double z = grid.cell_centroids[dim*c + dim - 1]; + below.reserve(number_of_cells); + above.reserve(number_of_cells); + for (int c = 0; c < number_of_cells; ++c) { + const double z = UgGridHelpers + ::getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, + dimension), + dimension-1); if (z > depth) { below.push_back(c); } else { @@ -76,8 +81,10 @@ namespace Opm // Initialize saturations so that there is water below woc, // and oil above. // If invert is true, water is instead above, oil below. - template - void initWaterOilContact(const UnstructuredGrid& grid, + template + void initWaterOilContact(int number_of_cells, + T begin_cell_centroids, + int dimensions, const Props& props, const double woc, const WaterInit waterinit, @@ -93,10 +100,12 @@ namespace Opm // } switch (waterinit) { case WaterBelow: - cellsBelowAbove(grid, woc, water, oil); + cellsBelowAbove(number_of_cells, begin_cell_centroids, dimensions, + woc, water, oil); break; case WaterAbove: - cellsBelowAbove(grid, woc, oil, water); + cellsBelowAbove(number_of_cells, begin_cell_centroids, dimensions, + woc, oil, water); } // Set saturations. state.setFirstSat(oil, props, State::MinSat); @@ -115,8 +124,10 @@ namespace Opm // Note that by manipulating the densities parameter, // it is possible to initialise with water on top instead of // on the bottom etc. - template - void initHydrostaticPressure(const UnstructuredGrid& grid, + template + void initHydrostaticPressure(int number_of_cells, + T begin_cell_centroids, + int dimensions, const double* densities, const double woc, const double gravity, @@ -125,14 +136,14 @@ namespace Opm State& state) { std::vector& p = state.pressure(); - const int num_cells = grid.number_of_cells; - const int dim = grid.dimensions; // Compute pressure at woc const double rho_datum = datum_z > woc ? densities[0] : densities[1]; const double woc_p = datum_p + (woc - datum_z)*gravity*rho_datum; - for (int c = 0; c < num_cells; ++c) { + for (int c = 0; c < number_of_cells; ++c) { // Compute pressure as delta from woc pressure. - const double z = grid.cell_centroids[dim*c + dim - 1]; + const double z = UgGridHelpers:: + getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions), + dimensions-1); const double rho = z > woc ? densities[0] : densities[1]; p[c] = woc_p + (z - woc)*gravity*rho; } @@ -141,8 +152,10 @@ namespace Opm // Facade to initHydrostaticPressure() taking a property object, // for similarity to initHydrostaticPressure() for compressible fluids. - template - void initHydrostaticPressure(const UnstructuredGrid& grid, + template + void initHydrostaticPressure(int number_of_cells, + T begin_cell_centroids, + int dimensions, const IncompPropertiesInterface& props, const double woc, const double gravity, @@ -151,7 +164,8 @@ namespace Opm State& state) { const double* densities = props.density(); - initHydrostaticPressure(grid, densities, woc, gravity, datum_z, datum_p, state); + initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions, + densities, woc, gravity, datum_z, datum_p, state); } @@ -183,8 +197,10 @@ namespace Opm // where rho is the oil density above the woc, water density // below woc. Note that even if there is (immobile) water in // the oil zone, it does not contribute to the pressure there. - template - void initHydrostaticPressure(const UnstructuredGrid& grid, + template + void initHydrostaticPressure(int number_of_cells, + T begin_cell_centroids, + int dimensions, const BlackoilPropertiesInterface& props, const double woc, const double gravity, @@ -195,11 +211,11 @@ namespace Opm assert(props.numPhases() == 2); // Obtain max and min z for which we will need to compute p. - const int num_cells = grid.number_of_cells; - const int dim = grid.dimensions; double zlim[2] = { 1e100, -1e100 }; - for (int c = 0; c < num_cells; ++c) { - const double z = grid.cell_centroids[dim*c + dim - 1]; + for (int c = 0; c < number_of_cells; ++c) { + const double z = UgGridHelpers:: + getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions), + dimensions-1);; zlim[0] = std::min(zlim[0], z); zlim[1] = std::max(zlim[1], z); } @@ -268,30 +284,42 @@ namespace Opm // Evaluate pressure at each cell centroid. std::vector& p = state.pressure(); - for (int c = 0; c < num_cells; ++c) { - const double z = grid.cell_centroids[dim*c + dim - 1]; + for (int c = 0; c < number_of_cells; ++c) { + const double z = UgGridHelpers:: + getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions), + dimensions-1); p[c] = press(z); } } // Initialize face pressures to distance-weighted average of adjacent cell pressures. - template - void initFacePressure(const UnstructuredGrid& grid, + template + void initFacePressure(int dimensions, + int number_of_faces, + FaceCells face_cells, + FCI begin_face_centroids, + CCI begin_cell_centroids, State& state) { - const int dim = grid.dimensions; const std::vector& cp = state.pressure(); std::vector& fp = state.facepressure(); - for (int f = 0; f < grid.number_of_faces; ++f) { + for (int f = 0; f < number_of_faces; ++f) { double dist[2] = { 0.0, 0.0 }; double press[2] = { 0.0, 0.0 }; int bdy_idx = -1; for (int j = 0; j < 2; ++j) { - const int c = grid.face_cells[2*f + j]; + const int c = face_cells(f, j); if (c >= 0) { dist[j] = 0.0; - for (int dd = 0; dd < dim; ++dd) { - double diff = grid.face_centroids[dim*f + dd] - grid.cell_centroids[dim*c + dd]; + for (int dd = 0; dd < dimensions; ++dd) { + double diff = UgGridHelpers:: + getCoordinate(UgGridHelpers::increment(begin_face_centroids, f, + dimensions), + dd) + - UgGridHelpers:: + getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, + dimensions), + dd); dist[j] += diff*diff; } dist[j] = std::sqrt(dist[j]); @@ -320,11 +348,32 @@ namespace Opm const double gravity, State& state) { + initStateBasic(grid.number_of_cells, grid.global_cell, grid.cartdims, + grid.number_of_faces, UgGridHelpers::faceCells(grid), + grid.face_centroids, grid.cell_centroids, + grid.dimensions, props, param, + gravity, state); + } + + template + void initStateBasic(int number_of_cells, + const int* global_cell, + const int* cartdims, + int number_of_faces, + FaceCells face_cells, + FCI begin_face_centroids, + CCI begin_cell_centroids, + int dimensions, + const IncompPropertiesInterface& props, + const parameter::ParameterGroup& param, + const double gravity, + State& state) +{ const int num_phases = props.numPhases(); if (num_phases != 2) { OPM_THROW(std::runtime_error, "initStateTwophaseBasic(): currently handling only two-phase scenarios."); } - state.init(grid, num_phases); + state.init(number_of_cells, number_of_faces, num_phases); const int num_cells = props.numCells(); // By default: initialise water saturation to minimum everywhere. std::vector all_cells(num_cells); @@ -338,11 +387,9 @@ namespace Opm // Initialise water saturation to max in the 'left' part. std::vector left_cells; left_cells.reserve(num_cells/2); - const int *glob_cell = grid.global_cell; - const int* cd = grid.cartdims; for (int cell = 0; cell < num_cells; ++cell) { - const int gc = glob_cell == 0 ? cell : glob_cell[cell]; - bool left = (gc % cd[0]) < cd[0]/2; + const int gc = global_cell == 0 ? cell : global_cell[cell]; + bool left = (gc % cartdims[0]) < cartdims[0]/2; if (left) { left_cells.push_back(cell); } @@ -355,30 +402,34 @@ namespace Opm if (gravity == 0.0) { std::cout << "**** Warning: running gravity segregation scenario, but gravity is zero." << std::endl; } - if (grid.cartdims[2] <= 1) { + if (cartdims[2] <= 1) { std::cout << "**** Warning: running gravity segregation scenario, which expects nz > 1." << std::endl; } // Initialise water saturation to max *above* water-oil contact. const double woc = param.get("water_oil_contact"); - initWaterOilContact(grid, props, woc, WaterAbove, state); + initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions, + props, woc, WaterAbove, state); // Initialise pressure to hydrostatic state. const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; double dens[2] = { props.density()[1], props.density()[0] }; - initHydrostaticPressure(grid, dens, woc, gravity, woc, ref_p, state); + initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions, + dens, woc, gravity, woc, ref_p, state); } else if (param.has("water_oil_contact")) { // Warn against error-prone usage. if (gravity == 0.0) { std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl; } - if (grid.cartdims[2] <= 1) { + if (cartdims[2] <= 1) { std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl; } // Initialise water saturation to max below water-oil contact. const double woc = param.get("water_oil_contact"); - initWaterOilContact(grid, props, woc, WaterBelow, state); + initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions, + props, woc, WaterBelow, state); // Initialise pressure to hydrostatic state. const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; - initHydrostaticPressure(grid, props.density(), woc, gravity, woc, ref_p, state); + initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions, + props.density(), woc, gravity, woc, ref_p, state); } else if (param.has("init_saturation")) { // Initialise water saturation to init_saturation parameter. const double init_saturation = param.get("init_saturation"); @@ -390,20 +441,23 @@ namespace Opm const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; const double rho = props.density()[0]*init_saturation + props.density()[1]*(1.0 - init_saturation); const double dens[2] = { rho, rho }; - const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1]; - initHydrostaticPressure(grid, dens, ref_z, gravity, ref_z, ref_p, state); + const double ref_z = UgGridHelpers::getCoordinate(begin_cell_centroids, dimensions - 1); + initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions, + dens, ref_z, gravity, ref_z, ref_p, state); } else { // Use default: water saturation is minimum everywhere. // Initialise pressure to hydrostatic state. const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; const double rho = props.density()[1]; const double dens[2] = { rho, rho }; - const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1]; - initHydrostaticPressure(grid, dens, ref_z, gravity, ref_z, ref_p, state); + const double ref_z = UgGridHelpers::getCoordinate(begin_cell_centroids, dimensions - 1); + initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions, + dens, ref_z, gravity, ref_z, ref_p, state); } // Finally, init face pressures. - initFacePressure(grid, state); + initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids, + begin_cell_centroids, state); } @@ -414,13 +468,33 @@ namespace Opm const parameter::ParameterGroup& param, const double gravity, State& state) + { + initStateBasic(grid.number_of_cells, grid.global_cell, grid.cartdims, + grid.number_of_faces, UgGridHelpers::faceCells(grid), + grid.face_centroids, grid.cell_centroids, grid.dimensions, + props, param, gravity, state); + } + + template + void initStateBasic(int number_of_cells, + const int* global_cell, + const int* cartdims, + int number_of_faces, + FaceCells face_cells, + FCI begin_face_centroids, + CCI begin_cell_centroids, + int dimensions, + const BlackoilPropertiesInterface& props, + const parameter::ParameterGroup& param, + const double gravity, + State& state) { // TODO: Refactor to exploit similarity with IncompProp* case. const int num_phases = props.numPhases(); if (num_phases != 2) { OPM_THROW(std::runtime_error, "initStateTwophaseBasic(): currently handling only two-phase scenarios."); } - state.init(grid, num_phases); + state.init(number_of_cells, number_of_faces, num_phases); const int num_cells = props.numCells(); // By default: initialise water saturation to minimum everywhere. std::vector all_cells(num_cells); @@ -433,11 +507,9 @@ namespace Opm // Initialise water saturation to max in the 'left' part. std::vector left_cells; left_cells.reserve(num_cells/2); - const int *glob_cell = grid.global_cell; - const int* cd = grid.cartdims; for (int cell = 0; cell < num_cells; ++cell) { - const int gc = glob_cell == 0 ? cell : glob_cell[cell]; - bool left = (gc % cd[0]) < cd[0]/2; + const int gc = global_cell == 0 ? cell : global_cell[cell]; + bool left = (gc % cartdims[0]) < cartdims[0]/2; if (left) { left_cells.push_back(cell); } @@ -450,26 +522,30 @@ namespace Opm if (gravity == 0.0) { std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl; } - if (grid.cartdims[2] <= 1) { + if (cartdims[2] <= 1) { std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl; } // Initialise water saturation to max below water-oil contact. const double woc = param.get("water_oil_contact"); - initWaterOilContact(grid, props, woc, WaterBelow, state); + initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions, + props, woc, WaterBelow, state); // Initialise pressure to hydrostatic state. const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; - initHydrostaticPressure(grid, props, woc, gravity, woc, ref_p, state); + initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions, + props, woc, gravity, woc, ref_p, state); } else { // Use default: water saturation is minimum everywhere. // Initialise pressure to hydrostatic state. const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; - const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1]; + const double ref_z = UgGridHelpers::getCoordinate(begin_cell_centroids, dimensions - 1); const double woc = -1e100; - initHydrostaticPressure(grid, props, woc, gravity, ref_z, ref_p, state); + initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions, + props, woc, gravity, ref_z, ref_p, state); } // Finally, init face pressures. - initFacePressure(grid, state); + initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids, + begin_cell_centroids, state); } @@ -577,6 +653,26 @@ namespace Opm const EclipseGridParser& deck, const double gravity, State& state) + { + initStateFromDeck(grid.number_of_cells, grid.global_cell, + grid.number_of_faces, UgGridHelpers::faceCells(grid), + grid.face_centroids, grid.cell_centroids, + grid.dimensions, props, deck, + gravity, state); + } + + template + void initStateFromDeck(int number_of_cells, + const int* global_cell, + int number_of_faces, + FaceCells face_cells, + FCI begin_face_centroids, + CCI begin_cell_centroids, + int dimensions, + const Props& props, + const EclipseGridParser& deck, + const double gravity, + State& state) { const int num_phases = props.numPhases(); const PhaseUsage pu = phaseUsageFromDeck(deck); @@ -584,7 +680,7 @@ namespace Opm OPM_THROW(std::runtime_error, "initStateFromDeck(): user specified property object with " << num_phases << " phases, " "found " << pu.num_phases << " phases in deck."); } - state.init(grid, num_phases); + state.init(number_of_cells, number_of_faces, num_phases); if (deck.hasField("EQUIL")) { if (num_phases != 2) { OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios."); @@ -598,17 +694,18 @@ namespace Opm OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet."); } const double woc = equil.equil[0].water_oil_contact_depth_; - initWaterOilContact(grid, props, woc, WaterBelow, state); + initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions, + props, woc, WaterBelow, state); // Set pressure depending on densities and depths. const double datum_z = equil.equil[0].datum_depth_; const double datum_p = equil.equil[0].datum_depth_pressure_; - initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state); + initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions, + props, woc, gravity, datum_z, datum_p, state); } else if (deck.hasField("PRESSURE")) { // Set saturations from SWAT/SGAS, pressure from PRESSURE. std::vector& s = state.saturation(); std::vector& p = state.pressure(); const std::vector& p_deck = deck.getFloatingPointValue("PRESSURE"); - const int num_cells = grid.number_of_cells; if (num_phases == 2) { if (!pu.phase_used[BlackoilPhases::Aqua]) { // oil-gas: we require SGAS @@ -618,8 +715,8 @@ namespace Opm const std::vector& sg_deck = deck.getFloatingPointValue("SGAS"); const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; const int opos = pu.phase_pos[BlackoilPhases::Liquid]; - for (int c = 0; c < num_cells; ++c) { - int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + for (int c = 0; c < number_of_cells; ++c) { + int c_deck = (global_cell == NULL) ? c : global_cell[c]; s[2*c + gpos] = sg_deck[c_deck]; s[2*c + opos] = 1.0 - sg_deck[c_deck]; p[c] = p_deck[c_deck]; @@ -632,8 +729,8 @@ namespace Opm const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; const int nwpos = (wpos + 1) % 2; - for (int c = 0; c < num_cells; ++c) { - int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + for (int c = 0; c < number_of_cells; ++c) { + int c_deck = (global_cell == NULL) ? c : global_cell[c]; s[2*c + wpos] = sw_deck[c_deck]; s[2*c + nwpos] = 1.0 - sw_deck[c_deck]; p[c] = p_deck[c_deck]; @@ -649,8 +746,8 @@ namespace Opm const int opos = pu.phase_pos[BlackoilPhases::Liquid]; const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); const std::vector& sg_deck = deck.getFloatingPointValue("SGAS"); - for (int c = 0; c < num_cells; ++c) { - int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + for (int c = 0; c < number_of_cells; ++c) { + int c_deck = (global_cell == NULL) ? c : global_cell[c]; s[3*c + wpos] = sw_deck[c_deck]; s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]); s[3*c + gpos] = sg_deck[c_deck]; @@ -664,7 +761,8 @@ namespace Opm } // Finally, init face pressures. - initFacePressure(grid, state); + initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids, + begin_cell_centroids, state); } /// Initialize surface volume from pressure and saturation by z = As. @@ -674,10 +772,18 @@ namespace Opm void initBlackoilSurfvol(const UnstructuredGrid& grid, const Props& props, State& state) + { + initBlackoilSurfvol(grid.number_of_cells, props, state); + } + + template + void initBlackoilSurfvol(int number_of_cells, + const Props& props, + State& state) { state.surfacevol() = state.saturation(); const int np = props.numPhases(); - const int nc = grid.number_of_cells; + const int nc = number_of_cells; std::vector allA(nc*np*np); std::vector allcells(nc); for (int c = 0; c < nc; ++c) { @@ -711,6 +817,15 @@ namespace Opm const Props& props, State& state) { + initBlackoilSurfvolUsingRSorRV(grid.number_of_cells, props, state); + } + + template + void initBlackoilSurfvolUsingRSorRV(int number_of_cells, + const Props& props, + State& state) + { + if (props.numPhases() != 3) { OPM_THROW(std::runtime_error, "initBlackoilSurfvol() is only supported in three-phase simulations."); } @@ -722,27 +837,26 @@ namespace Opm const PhaseUsage pu = props.phaseUsage(); const int np = props.numPhases(); - const int nc = grid.number_of_cells; - std::vector allA_a(nc*np*np); - std::vector allA_l(nc*np*np); - std::vector allA_v(nc*np*np); + std::vector allA_a(number_of_cells*np*np); + std::vector allA_l(number_of_cells*np*np); + std::vector allA_v(number_of_cells*np*np); - std::vector allcells(nc); - std::vector z_init(nc*np); + std::vector allcells(number_of_cells); + std::vector z_init(number_of_cells*np); std::fill(z_init.begin(),z_init.end(),0.0); - for (int c = 0; c < nc; ++c) { + for (int c = 0; c < number_of_cells; ++c) { allcells[c] = c; } - std::vector capPressures(nc*np); - props.capPress(nc,&state.saturation()[0],&allcells[0],&capPressures[0],NULL); + std::vector capPressures(number_of_cells*np); + props.capPress(number_of_cells,&state.saturation()[0],&allcells[0],&capPressures[0],NULL); - std::vector Pw(nc); - std::vector Pg(nc); + std::vector Pw(number_of_cells); + std::vector Pg(number_of_cells); - for (int c = 0; c < nc; ++c){ + for (int c = 0; c < number_of_cells; ++c){ Pw[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Aqua]; Pg[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Vapour]; } @@ -752,7 +866,7 @@ namespace Opm // Water phase if(pu.phase_used[BlackoilPhases::Aqua]) - for (int c = 0; c < nc ; ++c){ + for (int c = 0; c < number_of_cells ; ++c){ for (int p = 0; p < np ; ++p){ if (p == BlackoilPhases::Aqua) z_tmp = 1; @@ -762,11 +876,11 @@ namespace Opm z_init[c*np + p] = z_tmp; } } - props.matrix(nc, &state.pressure()[0], &z_init[0], &allcells[0], &allA_a[0], 0); + props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_a[0], 0); // Liquid phase if(pu.phase_used[BlackoilPhases::Liquid]){ - for (int c = 0; c < nc ; ++c){ + for (int c = 0; c < number_of_cells ; ++c){ for (int p = 0; p < np ; ++p){ if(p == BlackoilPhases::Vapour){ if(state.saturation()[np*c + p] > 0) @@ -784,10 +898,10 @@ namespace Opm } } } - props.matrix(nc, &state.pressure()[0], &z_init[0], &allcells[0], &allA_l[0], 0); + props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_l[0], 0); if(pu.phase_used[BlackoilPhases::Vapour]){ - for (int c = 0; c < nc ; ++c){ + for (int c = 0; c < number_of_cells ; ++c){ for (int p = 0; p < np ; ++p){ if(p == BlackoilPhases::Liquid){ if(state.saturation()[np*c + p] > 0) @@ -805,9 +919,9 @@ namespace Opm } } } - props.matrix(nc, &state.pressure()[0], &z_init[0], &allcells[0], &allA_v[0], 0); + props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_v[0], 0); - for (int c = 0; c < nc; ++c) { + for (int c = 0; c < number_of_cells; ++c) { // Using z = As double* z = &state.surfacevol()[c*np]; const double* s = &state.saturation()[c*np]; @@ -839,24 +953,43 @@ namespace Opm const double gravity, State& state) { - initStateFromDeck(grid, props, deck, gravity, state); + initBlackoilStateFromDeck(grid.number_of_cells, grid.global_cell, + grid.number_of_faces, UgGridHelpers::faceCells(grid), + grid.face_centroids, grid.cell_centroids,grid.dimensions, + props, deck, gravity, state); + } + + template + void initBlackoilStateFromDeck(int number_of_cells, + const int* global_cell, + int number_of_faces, + FaceCells face_cells, + FCI begin_face_centroids, + CCI begin_cell_centroids, + int dimensions, + const Props& props, + const EclipseGridParser& deck, + const double gravity, + State& state) + { + initStateFromDeck(number_of_cells, global_cell, number_of_faces, + face_cells, begin_face_centroids, begin_cell_centroids, + dimensions, props, deck, gravity, state); if (deck.hasField("RS")) { const std::vector& rs_deck = deck.getFloatingPointValue("RS"); - const int num_cells = grid.number_of_cells; - for (int c = 0; c < num_cells; ++c) { - int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + for (int c = 0; c < number_of_cells; ++c) { + int c_deck = (global_cell == NULL) ? c : global_cell[c]; state.gasoilratio()[c] = rs_deck[c_deck]; } - initBlackoilSurfvolUsingRSorRV(grid, props, state); + initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state); computeSaturation(props,state); } else if (deck.hasField("RV")){ const std::vector& rv_deck = deck.getFloatingPointValue("RV"); - const int num_cells = grid.number_of_cells; - for (int c = 0; c < num_cells; ++c) { - int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + for (int c = 0; c < number_of_cells; ++c) { + int c_deck = (global_cell == NULL) ? c : global_cell[c]; state.rv()[c] = rv_deck[c_deck]; } - initBlackoilSurfvolUsingRSorRV(grid, props, state); + initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state); computeSaturation(props,state); } diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index 92dffaef9..4b96ead72 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -45,15 +45,10 @@ namespace Opm const double* porosity, std::vector& porevol) { - int num_cells = grid.number_of_cells; - porevol.resize(num_cells); - std::transform(porosity, porosity + num_cells, - grid.cell_volumes, - porevol.begin(), - std::multiplies()); + computePorevolume(grid.number_of_cells, grid.cell_volumes, + porosity, porevol); } - /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects. /// @param[in] grid a grid /// @param[in] porosity array of grid.number_of_cells porosity values @@ -66,13 +61,11 @@ namespace Opm const std::vector& pressure, std::vector& porevol) { - int num_cells = grid.number_of_cells; - porevol.resize(num_cells); - for (int i = 0; i < num_cells; ++i) { - porevol[i] = porosity[i]*grid.cell_volumes[i]*rock_comp.poroMult(pressure[i]); - } + computePorevolume(grid.number_of_cells, grid.cell_volumes, porosity, rock_comp, pressure, + porevol); } + /// @brief Computes porosity of all cells in a grid, with rock compressibility effects. /// @param[in] grid a grid /// @param[in] porosity_standard array of grid.number_of_cells porosity values (at standard conditions) @@ -401,24 +394,16 @@ namespace Opm const std::vector& face_flux, std::vector& cell_velocity) { - const int dim = grid.dimensions; - cell_velocity.clear(); - cell_velocity.resize(grid.number_of_cells*dim, 0.0); - for (int face = 0; face < grid.number_of_faces; ++face) { - int c[2] = { grid.face_cells[2*face], grid.face_cells[2*face + 1] }; - const double* fc = &grid.face_centroids[face*dim]; - double flux = face_flux[face]; - for (int i = 0; i < 2; ++i) { - if (c[i] >= 0) { - const double* cc = &grid.cell_centroids[c[i]*dim]; - for (int d = 0; d < dim; ++d) { - double v_contrib = fc[d] - cc[d]; - v_contrib *= flux/grid.cell_volumes[c[i]]; - cell_velocity[c[i]*dim + d] += (i == 0) ? v_contrib : -v_contrib; - } - } - } - } + estimateCellVelocity(grid.number_of_cells, + grid.number_of_faces, + grid.face_centroids, + UgGridHelpers::faceCells(grid), + grid.cell_centroids, + grid.cell_volumes, + grid.dimensions, + face_flux, + cell_velocity); + } /// Extract a vector of water saturations from a vector of @@ -490,44 +475,8 @@ namespace Opm const double* densities, const double gravity, const bool per_grid_cell, std::vector& wdp) { - const int nw = wells.number_of_wells; - const size_t np = per_grid_cell ? - saturations.size()/grid.number_of_cells - : saturations.size()/wells.well_connpos[nw]; - // Simple for now: - for (int i = 0; i < nw; i++) { - double depth_ref = wells.depth_ref[i]; - for (int j = wells.well_connpos[i]; j < wells.well_connpos[i + 1]; j++) { - int cell = wells.well_cells[j]; - - // Is this correct wrt. depth_ref? - double cell_depth = grid.cell_centroids[3 * cell + 2]; - - double saturation_sum = 0.0; - for (size_t p = 0; p < np; p++) { - if (!per_grid_cell) { - saturation_sum += saturations[j * np + p]; - } else { - saturation_sum += saturations[np * cell + p]; - } - } - if (saturation_sum == 0) { - saturation_sum = 1.0; - } - double density = 0.0; - for (size_t p = 0; p < np; p++) { - if (!per_grid_cell) { - density += saturations[j * np + p] * densities[p] / saturation_sum; - } else { - // Is this a smart way of doing it? - density += saturations[np * cell + p] * densities[p] / saturation_sum; - } - } - - // Is the sign correct? - wdp.push_back(density * (cell_depth - depth_ref) * gravity); - } - } + computeWDP(wells, grid.number_of_cells, grid.cell_centroids, saturations, densities, + gravity, per_grid_cell, wdp); } diff --git a/opm/core/utility/miscUtilities.hpp b/opm/core/utility/miscUtilities.hpp index 05bfaf2ac..1362e84d7 100644 --- a/opm/core/utility/miscUtilities.hpp +++ b/opm/core/utility/miscUtilities.hpp @@ -41,6 +41,16 @@ namespace Opm const double* porosity, std::vector& porevol); + /// @brief Computes pore volume of all cells in a grid. + /// @param[in] number_of_cells The number of cells of the grid. + /// @param[in] begin_cell_volume Iterator to the volume of the first cell. + /// @param[in] porosity array of grid.number_of_cells porosity values + /// @param[out] porevol the pore volume by cell. + template + void computePorevolume(int number_of_cells, + T begin_cell_volume, + const double* porosity, + std::vector& porevol); /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects. /// @param[in] grid a grid @@ -54,6 +64,21 @@ namespace Opm const std::vector& pressure, std::vector& porevol); + /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects. + /// @param[in] number_of_cells The number of cells of the grid. + /// @param[in] Pointer to/ Iterator at the first cell volume. + /// @param[in] porosity array of grid.number_of_cells porosity values + /// @param[in] rock_comp rock compressibility properties + /// @param[in] pressure pressure by cell + /// @param[out] porevol the pore volume by cell. + template + void computePorevolume(int number_of_cells, + T begin_cell_volume, + const double* porosity, + const RockCompressibility& rock_comp, + const std::vector& pressure, + std::vector& porevol); + /// @brief Computes porosity of all cells in a grid, with rock compressibility effects. /// @param[in] grid a grid /// @param[in] porosity_standard array of grid.number_of_cells porosity values (at reference presure) @@ -188,6 +213,26 @@ namespace Opm const std::vector& face_flux, std::vector& cell_velocity); + /// @brief Estimates a scalar cell velocity from face fluxes. + /// @param[in] number_of_cells The number of cells of the grid + /// @param[in] number_of_faces The number of cells of the grid + /// @param[in] begin_face_centroids Iterator pointing to first face centroid. + /// @param[in] face_cells Mapping from faces to connected cells. + /// @param[in] dimensions The dimensions of the grid. + /// @param[in] begin_cell_centroids Iterator pointing to first cell centroid. + /// @param[in] face_flux signed per-face fluxes + /// @param[out] cell_velocity the estimated velocities. + template + void estimateCellVelocity(int number_of_cells, + int number_of_faces, + FC begin_face_centroids, + FC1 face_cells, + CC begin_cell_centroids, + CV begin_cell_volumes, + int dimension, + const std::vector& face_flux, + std::vector& cell_velocity); + /// Extract a vector of water saturations from a vector of /// interleaved water and oil saturations. void toWaterSat(const std::vector& sboth, @@ -218,6 +263,24 @@ namespace Opm const double* densities, const double gravity, const bool per_grid_cell, std::vector& wdp); + /// Computes the WDP for each well. + /// \param[in] wells Wells that need their wdp calculated. + /// \param[in] number_of_cells The number of cells in the grid. + /// \param[in] begin_cell_centroids Pointer/Iterator to the first cell centroid. + /// \param[in] saturations A vector of weights for each cell for each phase + /// in the grid (or well, see per_grid_cell parameter). So for cell i, + /// saturations[i*densities.size() + p] should give the weight + /// of phase p in cell i. + /// \param[in] densities Density for each phase. + /// \param[out] wdp Will contain, for each well, the wdp of the well. + /// \param[in] per_grid_cell Whether or not the saturations are per grid cell or per + /// well cell. + template + void computeWDP(const Wells& wells, int number_of_cells, T begin_cell_centroids, + const std::vector& saturations, + const double* densities, const double gravity, const bool per_grid_cell, + std::vector& wdp); + /// Computes (sums) the flow rate for each well. /// \param[in] wells The wells for which the flow rate should be computed. /// \param[in] flow_rates_per_cell Flow rates per well cells. Should ordered the same way as @@ -320,4 +383,5 @@ namespace Opm } // namespace Opm +#include "miscUtilities_impl.hpp" #endif // OPM_MISCUTILITIES_HEADER_INCLUDED diff --git a/opm/core/utility/miscUtilities_impl.hpp b/opm/core/utility/miscUtilities_impl.hpp new file mode 100644 index 000000000..25bc56111 --- /dev/null +++ b/opm/core/utility/miscUtilities_impl.hpp @@ -0,0 +1,123 @@ +#include +#include +#include + +namespace Opm +{ + /// @brief Estimates a scalar cell velocity from face fluxes. + /// @param[in] number_of_cells The number of cells of the grid + /// @param[in] begin_face_centroids Iterator pointing to first face centroid. + /// @param[in] face_cells Mapping from faces to connected cells. + /// @param[in] dimensions The dimensions of the grid. + /// @param[in] begin_cell_centroids Iterator pointing to first cell centroid. + /// @param[in] face_flux signed per-face fluxes + /// @param[out] cell_velocity the estimated velocities. + template + void estimateCellVelocity(int number_of_cells, + int number_of_faces, + FC begin_face_centroids, + FC1 face_cells, + CC begin_cell_centroids, + CV begin_cell_volumes, + int dimension, + const std::vector& face_flux, + std::vector& cell_velocity) + { + cell_velocity.clear(); + cell_velocity.resize(number_of_cells*dimension, 0.0); + for (int face = 0; face < number_of_faces; ++face) { + int c[2] = { face_cells(face, 0), face_cells(face, 1) }; + FC fc = UgGridHelpers::increment(begin_face_centroids, face, dimension); + double flux = face_flux[face]; + for (int i = 0; i < 2; ++i) { + if (c[i] >= 0) { + CC cc = UgGridHelpers::increment(begin_cell_centroids, c[i], dimension); + for (int d = 0; d < dimension; ++d) { + double v_contrib = UgGridHelpers::getCoordinate(fc, d) - UgGridHelpers::getCoordinate(cc, d); + v_contrib *= flux/begin_cell_volumes[c[i]]; + cell_velocity[c[i]*dimension + d] += (i == 0) ? v_contrib : -v_contrib; + } + } + } + } + } + + template + void computePorevolume(int number_of_cells, + T begin_cell_volume, + const double* porosity, + std::vector& porevol) + { + porevol.resize(number_of_cells); + std::transform(porosity, porosity + number_of_cells, + begin_cell_volume, + porevol.begin(), + std::multiplies()); + } + + /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects. + /// @param[in] number_of_cells The number of cells of the grid. + /// @param[in] porosity array of grid.number_of_cells porosity values + /// @param[in] rock_comp rock compressibility properties + /// @param[in] pressure pressure by cell + /// @param[out] porevol the pore volume by cell. + template + void computePorevolume(int number_of_cells, + T begin_cell_volumes, + const double* porosity, + const RockCompressibility& rock_comp, + const std::vector& pressure, + std::vector& porevol) + { + porevol.resize(number_of_cells); + for (int i = 0; i < number_of_cells; ++i) { + porevol[i] = porosity[i]*begin_cell_volumes[i]*rock_comp.poroMult(pressure[i]); + } + } + + template + void computeWDP(const Wells& wells, int number_of_cells, T begin_cell_centroids, const std::vector& saturations, + const double* densities, const double gravity, const bool per_grid_cell, + std::vector& wdp) + { + const int nw = wells.number_of_wells; + const size_t np = per_grid_cell ? + saturations.size()/number_of_cells + : saturations.size()/wells.well_connpos[nw]; + // Simple for now: + for (int i = 0; i < nw; i++) { + double depth_ref = wells.depth_ref[i]; + for (int j = wells.well_connpos[i]; j < wells.well_connpos[i + 1]; j++) { + int cell = wells.well_cells[j]; + + // Is this correct wrt. depth_ref? + double cell_depth = UgGridHelpers + ::getCoordinate(UgGridHelpers::increment(begin_cell_centroids, cell, 3), 2); + + double saturation_sum = 0.0; + for (size_t p = 0; p < np; p++) { + if (!per_grid_cell) { + saturation_sum += saturations[j * np + p]; + } else { + saturation_sum += saturations[np * cell + p]; + } + } + if (saturation_sum == 0) { + saturation_sum = 1.0; + } + double density = 0.0; + for (size_t p = 0; p < np; p++) { + if (!per_grid_cell) { + density += saturations[j * np + p] * densities[p] / saturation_sum; + } else { + // Is this a smart way of doing it? + density += saturations[np * cell + p] * densities[p] / saturation_sum; + } + } + + // Is the sign correct? + wdp.push_back(density * (cell_depth - depth_ref) * gravity); + } + } + } +} diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 7df448cf8..24ad6c288 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -26,13 +26,11 @@ #include #include #include -#include #include #include #include -#include #include #include #include @@ -44,17 +42,13 @@ // Helper structs and functions for the implementation. -namespace +namespace WellsManagerDetail { namespace ProductionControl { - enum Mode { ORAT, WRAT, GRAT, - LRAT, CRAT, RESV, - BHP , THP , GRUP }; - namespace Details { std::map init_mode_map() { @@ -122,8 +116,6 @@ namespace namespace InjectionControl { - enum Mode { RATE, RESV, BHP, - THP, GRUP }; namespace Details { std::map @@ -177,28 +169,6 @@ namespace } // namespace InjectionControl - - std::array getCubeDim(const UnstructuredGrid& grid, int cell) - { - using namespace std; - std::array cube; - int num_local_faces = grid.cell_facepos[cell + 1] - grid.cell_facepos[cell]; - vector x(num_local_faces); - vector y(num_local_faces); - vector z(num_local_faces); - for (int lf=0; lf cartesian_to_compressed; - setupCompressedToCartesian(grid, cartesian_to_compressed); + setupCompressedToCartesian(UgGridHelpers::globalCell(grid), + UgGridHelpers::numCells(grid), cartesian_to_compressed); // Obtain phase usage data. PhaseUsage pu = phaseUsageFromDeck(eclipseState); @@ -305,7 +277,12 @@ namespace Opm well_names.reserve(wells.size()); well_data.reserve(wells.size()); - createWellsFromSpecs(wells, timeStep, grid, well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability); + createWellsFromSpecs(wells, timeStep, UgGridHelpers::cell2Faces(grid), + UgGridHelpers::cartDims(grid), + UgGridHelpers::beginFaceCentroids(grid), + UgGridHelpers::beginCellCentroids(grid), + UgGridHelpers::dimensions(grid), + well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability); setupWellControls(wells, timeStep, well_names, pu); { @@ -352,586 +329,16 @@ namespace Opm /// Construct wells from deck. WellsManager::WellsManager(const Opm::EclipseGridParser& deck, const UnstructuredGrid& grid, - const double* permeability) - : w_(0) + const double* permeability, + bool checkCellExistence) + : w_(0), checkCellExistence_(checkCellExistence) { - if (grid.dimensions != 3) { - OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional."); - } - // NOTE: Implementation copied and modified from dune-porsol's class BlackoilWells. - std::vector keywords; - keywords.push_back("WELSPECS"); - keywords.push_back("COMPDAT"); -// keywords.push_back("WELTARG"); - if (!deck.hasFields(keywords)) { - OPM_MESSAGE("Missing well keywords in deck, initializing no wells."); - return; - } - if (!(deck.hasField("WCONINJE") || deck.hasField("WCONPROD")) ) { - OPM_THROW(std::runtime_error, "Needed field is missing in file"); - } - - // Obtain phase usage data. - PhaseUsage pu = phaseUsageFromDeck(deck); - - // These data structures will be filled in this constructor, - // then used to initialize the Wells struct. - std::vector well_names; - std::vector well_data; - std::vector > wellperf_data; - - // For easy lookup: - std::map well_names_to_index; - typedef std::map::const_iterator WNameIt; - - // Get WELSPECS data. - // It is allowed to have multiple lines corresponding to - // the same well, in which case the last one encountered - // is the one used. - const WELSPECS& welspecs = deck.getWELSPECS(); - const int num_welspecs = welspecs.welspecs.size(); - well_names.reserve(num_welspecs); - well_data.reserve(num_welspecs); - for (int w = 0; w < num_welspecs; ++w) { - // First check if this well has already been encountered. - // If so, we modify it's data instead of appending a new well - // to the well_data and well_names vectors. - const std::string& name = welspecs.welspecs[w].name_; - const double refdepth = welspecs.welspecs[w].datum_depth_BHP_; - WNameIt wit = well_names_to_index.find(name); - if (wit == well_names_to_index.end()) { - // New well, append data. - well_names_to_index[welspecs.welspecs[w].name_] = well_data.size(); - well_names.push_back(name); - WellData wd; - // If negative (defaulted), set refdepth to a marker - // value, will be changed after getting perforation - // data to the centroid of the cell of the top well - // perforation. - wd.reference_bhp_depth = (refdepth < 0.0) ? -1e100 : refdepth; - wd.welspecsline = w; - well_data.push_back(wd); - } else { - // Existing well, change data. - const int wix = wit->second; - well_data[wix].reference_bhp_depth = (refdepth < 0.0) ? -1e100 : refdepth; - well_data[wix].welspecsline = w; - } - } - const int num_wells = well_data.size(); - wellperf_data.resize(num_wells); - - - // global_cell is a map from compressed cells to Cartesian grid cells. - // We must make the inverse lookup. - const int* global_cell = grid.global_cell; - const int* cpgdim = grid.cartdims; - std::map cartesian_to_compressed; - - if (global_cell) { - for (int i = 0; i < grid.number_of_cells; ++i) { - cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); - } - } - else { - for (int i = 0; i < grid.number_of_cells; ++i) { - cartesian_to_compressed.insert(std::make_pair(i, i)); - } - } - - // Get COMPDAT data - // It is *not* allowed to have multiple lines corresponding to - // the same perforation! - const COMPDAT& compdat = deck.getCOMPDAT(); - const int num_compdat = compdat.compdat.size(); - for (int kw = 0; kw < num_compdat; ++kw) { - // Extract well name, or the part of the well name that - // comes before the '*'. - std::string name = compdat.compdat[kw].well_; - std::string::size_type len = name.find('*'); - if (len != std::string::npos) { - name = name.substr(0, len); - } - // Look for well with matching name. - bool found = false; - for (int wix = 0; wix < num_wells; ++wix) { - if (well_names[wix].compare(0,len, name) == 0) { // equal - // Extract corresponding WELSPECS defintion for - // purpose of default location specification. - const WelspecsLine& wspec = welspecs.welspecs[well_data[wix].welspecsline]; - - // We have a matching name. - int ix = compdat.compdat[kw].grid_ind_[0] - 1; - int jy = compdat.compdat[kw].grid_ind_[1] - 1; - int kz1 = compdat.compdat[kw].grid_ind_[2] - 1; - int kz2 = compdat.compdat[kw].grid_ind_[3] - 1; - - if (ix < 0) { - // Defaulted I location. Extract from WELSPECS. - ix = wspec.I_ - 1; - } - if (jy < 0) { - // Defaulted J location. Extract from WELSPECS. - jy = wspec.J_ - 1; - } - if (kz1 < 0) { - // Defaulted KZ1. Use top layer. - kz1 = 0; - } - if (kz2 < 0) { - // Defaulted KZ2. Use bottom layer. - kz2 = cpgdim[2] - 1; - } - - for (int kz = kz1; kz <= kz2; ++kz) { - int cart_grid_indx = ix + cpgdim[0]*(jy + cpgdim[1]*kz); - std::map::const_iterator cgit = - cartesian_to_compressed.find(cart_grid_indx); - if (cgit == cartesian_to_compressed.end()) { - OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << ix << ' ' << jy << ' ' - << kz << " not found in grid (well = " << name << ')'); - } - int cell = cgit->second; - PerfData pd; - pd.cell = cell; - if (compdat.compdat[kw].connect_trans_fac_ > 0.0) { - pd.well_index = compdat.compdat[kw].connect_trans_fac_; - } else { - double radius = 0.5*compdat.compdat[kw].diameter_; - if (radius <= 0.0) { - radius = 0.5*unit::feet; - OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); - } - std::array cubical = getCubeDim(grid, cell); - const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell]; - pd.well_index = computeWellIndex(radius, cubical, cell_perm, - compdat.compdat[kw].skin_factor_); - } - wellperf_data[wix].push_back(pd); - } - found = true; - break; - } - } - if (!found) { - OPM_THROW(std::runtime_error, "Undefined well name: " << compdat.compdat[kw].well_ - << " in COMPDAT"); - } - } - - // Set up reference depths that were defaulted. Count perfs. - int num_perfs = 0; - assert(grid.dimensions == 3); - for (int w = 0; w < num_wells; ++w) { - num_perfs += wellperf_data[w].size(); - if (well_data[w].reference_bhp_depth < 0.0) { - // It was defaulted. Set reference depth to minimum perforation depth. - double min_depth = 1e100; - int num_wperfs = wellperf_data[w].size(); - for (int perf = 0; perf < num_wperfs; ++perf) { - double depth = grid.cell_centroids[3*wellperf_data[w][perf].cell + 2]; - min_depth = std::min(min_depth, depth); - } - well_data[w].reference_bhp_depth = min_depth; - } - } - - // Create the well data structures. - w_ = create_wells(pu.num_phases, num_wells, num_perfs); - if (!w_) { - OPM_THROW(std::runtime_error, "Failed creating Wells struct."); - } - - // Classify wells - if (deck.hasField("WCONINJE")) { - const std::vector& lines = deck.getWCONINJE().wconinje; - for (size_t i = 0 ; i < lines.size(); ++i) { - const std::map::const_iterator it = well_names_to_index.find(lines[i].well_); - if (it != well_names_to_index.end()) { - const int well_index = it->second; - well_data[well_index].type = INJECTOR; - } else { - OPM_THROW(std::runtime_error, "Unseen well name: " << lines[i].well_ << " first seen in WCONINJE"); - } - } - } - if (deck.hasField("WCONPROD")) { - const std::vector& lines = deck.getWCONPROD().wconprod; - for (size_t i = 0; i < lines.size(); ++i) { - const std::map::const_iterator it = well_names_to_index.find(lines[i].well_); - if (it != well_names_to_index.end()) { - const int well_index = it->second; - well_data[well_index].type = PRODUCER; - } else { - OPM_THROW(std::runtime_error, "Unseen well name: " << lines[i].well_ << " first seen in WCONPROD"); - } - - } - } - - // Add wells. - for (int w = 0; w < num_wells; ++w) { - const int w_num_perf = wellperf_data[w].size(); - std::vector perf_cells(w_num_perf); - std::vector perf_prodind(w_num_perf); - for (int perf = 0; perf < w_num_perf; ++perf) { - perf_cells[perf] = wellperf_data[w][perf].cell; - perf_prodind[perf] = wellperf_data[w][perf].well_index; - } - const double* comp_frac = NULL; - // We initialize all wells with a null component fraction, - // and must (for injection wells) overwrite it later. - int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, w_num_perf, - comp_frac, &perf_cells[0], &perf_prodind[0], well_names[w].c_str(), w_); - if (!ok) { - OPM_THROW(std::runtime_error, "Failed adding well " << well_names[w] << " to Wells data structure."); - } - } - - // Get WCONINJE data, add injection controls to wells. - // It is allowed to have multiple lines corresponding to - // the same well, in which case the last one encountered - // is the one used. - if (deck.hasField("WCONINJE")) { - const WCONINJE& wconinjes = deck.getWCONINJE(); - const int num_wconinjes = wconinjes.wconinje.size(); - for (int kw = 0; kw < num_wconinjes; ++kw) { - const WconinjeLine& wci_line = wconinjes.wconinje[kw]; - // Extract well name, or the part of the well name that - // comes before the '*'. - std::string name = wci_line.well_; - std::string::size_type len = name.find('*'); - if (len != std::string::npos) { - name = name.substr(0, len); - } - bool well_found = false; - for (int wix = 0; wix < num_wells; ++wix) { - if (well_names[wix].compare(0,len, name) == 0) { //equal - well_found = true; - assert(well_data[wix].type == w_->type[wix]); - if (well_data[wix].type != INJECTOR) { - OPM_THROW(std::runtime_error, "Found WCONINJE entry for a non-injector well: " << well_names[wix]); - } - - // Add all controls that are present in well. - // First we must clear existing controls, in case the - // current WCONINJE line is modifying earlier controls. - clear_well_controls(wix, w_); - int ok = 1; - int control_pos[5] = { -1, -1, -1, -1, -1 }; - if (ok && wci_line.surface_flow_max_rate_ >= 0.0) { - control_pos[InjectionControl::RATE] = well_controls_get_num(w_->ctrls[wix]); - double distr[3] = { 0.0, 0.0, 0.0 }; - if (wci_line.injector_type_ == "WATER") { - distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; - } else if (wci_line.injector_type_ == "OIL") { - distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; - } else if (wci_line.injector_type_ == "GAS") { - distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; - } else { - OPM_THROW(std::runtime_error, "Injector type " << wci_line.injector_type_ << " not supported." - "WellsManager only supports WATER, OIL and GAS injector types."); - } - ok = append_well_controls(SURFACE_RATE, wci_line.surface_flow_max_rate_, - distr, wix, w_); - } - if (ok && wci_line.reservoir_flow_max_rate_ >= 0.0) { - control_pos[InjectionControl::RESV] = well_controls_get_num(w_->ctrls[wix]); - double distr[3] = { 0.0, 0.0, 0.0 }; - if (wci_line.injector_type_ == "WATER") { - distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; - } else if (wci_line.injector_type_ == "OIL") { - distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; - } else if (wci_line.injector_type_ == "GAS") { - distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; - } else { - OPM_THROW(std::runtime_error, "Injector type " << wci_line.injector_type_ << " not supported." - "WellsManager only supports WATER, OIL and GAS injector types."); - } - ok = append_well_controls(RESERVOIR_RATE, wci_line.reservoir_flow_max_rate_, - distr, wix, w_); - } - if (ok && wci_line.BHP_limit_ > 0.0) { - control_pos[InjectionControl::BHP] = well_controls_get_num(w_->ctrls[wix]); - ok = append_well_controls(BHP, wci_line.BHP_limit_, - NULL, wix, w_); - } - if (ok && wci_line.THP_limit_ > 0.0) { - OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[wix]); - } - if (!ok) { - OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[wix]); - } - InjectionControl::Mode mode = InjectionControl::mode(wci_line.control_mode_); - int cpos = control_pos[mode]; - if (cpos == -1 && mode != InjectionControl::GRUP) { - OPM_THROW(std::runtime_error, "Control for " << wci_line.control_mode_ << " not specified in well " << well_names[wix]); - } - // We need to check if the well is shut or not - if (wci_line.open_shut_flag_ == "SHUT") { - cpos = ~cpos; - } - set_current_control(wix, cpos, w_); - - // Set well component fraction. - double cf[3] = { 0.0, 0.0, 0.0 }; - if (wci_line.injector_type_[0] == 'W') { - if (!pu.phase_used[BlackoilPhases::Aqua]) { - OPM_THROW(std::runtime_error, "Water phase not used, yet found water-injecting well."); - } - cf[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; - } else if (wci_line.injector_type_[0] == 'O') { - if (!pu.phase_used[BlackoilPhases::Liquid]) { - OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-injecting well."); - } - cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; - } else if (wci_line.injector_type_[0] == 'G') { - if (!pu.phase_used[BlackoilPhases::Vapour]) { - OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-injecting well."); - } - cf[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; - } - std::copy(cf, cf + pu.num_phases, w_->comp_frac + wix*pu.num_phases); - } - } - if (!well_found) { - OPM_THROW(std::runtime_error, "Undefined well name: " << wci_line.well_ - << " in WCONINJE"); - } - } - } - - // Get WCONPROD data - // It is allowed to have multiple lines corresponding to - // the same well, in which case the last one encountered - // is the one used. - if (deck.hasField("WCONPROD")) { - const WCONPROD& wconprods = deck.getWCONPROD(); - const int num_wconprods = wconprods.wconprod.size(); - for (int kw = 0; kw < num_wconprods; ++kw) { - const WconprodLine& wcp_line = wconprods.wconprod[kw]; - std::string name = wcp_line.well_; - std::string::size_type len = name.find('*'); - if (len != std::string::npos) { - name = name.substr(0, len); - } - bool well_found = false; - for (int wix = 0; wix < num_wells; ++wix) { - if (well_names[wix].compare(0,len, name) == 0) { //equal - well_found = true; - assert(well_data[wix].type == w_->type[wix]); - if (well_data[wix].type != PRODUCER) { - OPM_THROW(std::runtime_error, "Found WCONPROD entry for a non-producer well: " << well_names[wix]); - } - // Add all controls that are present in well. - // First we must clear existing controls, in case the - // current WCONPROD line is modifying earlier controls. - clear_well_controls(wix, w_); - int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 }; - int ok = 1; - if (ok && wcp_line.oil_max_rate_ >= 0.0) { - if (!pu.phase_used[BlackoilPhases::Liquid]) { - OPM_THROW(std::runtime_error, "Oil phase not active and ORAT control specified."); - } - control_pos[ProductionControl::ORAT] = well_controls_get_num(w_->ctrls[wix]); - double distr[3] = { 0.0, 0.0, 0.0 }; - distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; - ok = append_well_controls(SURFACE_RATE, -wcp_line.oil_max_rate_, - distr, wix, w_); - } - if (ok && wcp_line.water_max_rate_ >= 0.0) { - if (!pu.phase_used[BlackoilPhases::Aqua]) { - OPM_THROW(std::runtime_error, "Water phase not active and WRAT control specified."); - } - control_pos[ProductionControl::WRAT] = well_controls_get_num(w_->ctrls[wix]); - double distr[3] = { 0.0, 0.0, 0.0 }; - distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; - ok = append_well_controls(SURFACE_RATE, -wcp_line.water_max_rate_, - distr, wix, w_); - } - if (ok && wcp_line.gas_max_rate_ >= 0.0) { - if (!pu.phase_used[BlackoilPhases::Vapour]) { - OPM_THROW(std::runtime_error, "Gas phase not active and GRAT control specified."); - } - control_pos[ProductionControl::GRAT] = well_controls_get_num(w_->ctrls[wix]); - double distr[3] = { 0.0, 0.0, 0.0 }; - distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; - ok = append_well_controls(SURFACE_RATE, -wcp_line.gas_max_rate_, - distr, wix, w_); - } - if (ok && wcp_line.liquid_max_rate_ >= 0.0) { - if (!pu.phase_used[BlackoilPhases::Aqua]) { - OPM_THROW(std::runtime_error, "Water phase not active and LRAT control specified."); - } - if (!pu.phase_used[BlackoilPhases::Liquid]) { - OPM_THROW(std::runtime_error, "Oil phase not active and LRAT control specified."); - } - control_pos[ProductionControl::LRAT] = well_controls_get_num(w_->ctrls[wix]); - double distr[3] = { 0.0, 0.0, 0.0 }; - distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; - distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; - ok = append_well_controls(SURFACE_RATE, -wcp_line.liquid_max_rate_, - distr, wix, w_); - } - if (ok && wcp_line.reservoir_flow_max_rate_ >= 0.0) { - control_pos[ProductionControl::RESV] = well_controls_get_num(w_->ctrls[wix]); - double distr[3] = { 1.0, 1.0, 1.0 }; - ok = append_well_controls(RESERVOIR_RATE, -wcp_line.reservoir_flow_max_rate_, - distr, wix, w_); - } - if (ok && wcp_line.BHP_limit_ > 0.0) { - control_pos[ProductionControl::BHP] = well_controls_get_num(w_->ctrls[wix]); - ok = append_well_controls(BHP, wcp_line.BHP_limit_, - NULL, wix, w_); - } - if (ok && wcp_line.THP_limit_ > 0.0) { - OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[wix]); - } - if (!ok) { - OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[wix]); - } - ProductionControl::Mode mode = ProductionControl::mode(wcp_line.control_mode_); - int cpos = control_pos[mode]; - if (cpos == -1 && mode != ProductionControl::GRUP) { - OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[wix]); - } - // If it's shut, we complement the cpos - if (wcp_line.open_shut_flag_ == "SHUT") { - cpos = ~cpos; // So we can easily retrieve the cpos later - } - set_current_control(wix, cpos, w_); - } - } - if (!well_found) { - OPM_THROW(std::runtime_error, "Undefined well name: " << wcp_line.well_ - << " in WCONPROD"); - } - } - } - - // Get WELTARG data - if (deck.hasField("WELTARG")) { - OPM_THROW(std::runtime_error, "We currently do not handle WELTARG."); - /* - const WELTARG& weltargs = deck.getWELTARG(); - const int num_weltargs = weltargs.weltarg.size(); - for (int kw = 0; kw < num_weltargs; ++kw) { - std::string name = weltargs.weltarg[kw].well_; - std::string::size_type len = name.find('*'); - if (len != std::string::npos) { - name = name.substr(0, len); - } - bool well_found = false; - for (int wix = 0; wix < num_wells; ++wix) { - if (well_names[wix].compare(0,len, name) == 0) { //equal - well_found = true; - well_data[wix].target = weltargs.weltarg[kw].new_value_; - break; - } - } - if (!well_found) { - OPM_THROW(std::runtime_error, "Undefined well name: " << weltargs.weltarg[kw].well_ - << " in WELTARG"); - } - } - */ - } - - // Debug output. -#define EXTRA_OUTPUT -#ifdef EXTRA_OUTPUT - /* - std::cout << "\t WELL DATA" << std::endl; - for(int i = 0; i< num_wells; ++i) { - std::cout << i << ": " << well_data[i].type << " " - << well_data[i].control << " " << well_data[i].target - << std::endl; - } - - std::cout << "\n\t PERF DATA" << std::endl; - for(int i=0; i< int(wellperf_data.size()); ++i) { - for(int j=0; j< int(wellperf_data[i].size()); ++j) { - std::cout << i << ": " << wellperf_data[i][j].cell << " " - << wellperf_data[i][j].well_index << std::endl; - } - } - */ -#endif - - if (deck.hasField("WELOPEN")) { - const WELOPEN& welopen = deck.getWELOPEN(); - for (size_t i = 0; i < welopen.welopen.size(); ++i) { - WelopenLine line = welopen.welopen[i]; - std::string wellname = line.well_; - std::map::const_iterator it = well_names_to_index.find(wellname); - if (it == well_names_to_index.end()) { - OPM_THROW(std::runtime_error, "Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS."); - } - const int index = it->second; - if (line.openshutflag_ == "SHUT") { - well_controls_shut_well( w_->ctrls[index] ); - } else if (line.openshutflag_ == "OPEN") { - well_controls_open_well( w_->ctrls[index] ); - } else { - OPM_THROW(std::runtime_error, "Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT."); - } - } - } - - // Build the well_collection_ well group hierarchy. - if (deck.hasField("GRUPTREE")) { - std::cout << "Found gruptree" << std::endl; - const GRUPTREE& gruptree = deck.getGRUPTREE(); - std::map::const_iterator it = gruptree.tree.begin(); - for( ; it != gruptree.tree.end(); ++it) { - well_collection_.addChild(it->first, it->second, deck); - } - } - for (size_t i = 0; i < welspecs.welspecs.size(); ++i) { - WelspecsLine line = welspecs.welspecs[i]; - well_collection_.addChild(line.name_, line.group_, deck); - } - - - - // Set the guide rates: - if (deck.hasField("WGRUPCON")) { - std::cout << "Found Wgrupcon" << std::endl; - WGRUPCON wgrupcon = deck.getWGRUPCON(); - const std::vector& lines = wgrupcon.wgrupcon; - std::cout << well_collection_.getLeafNodes().size() << std::endl; - for (size_t i = 0; i < lines.size(); i++) { - std::string name = lines[i].well_; - const int wix = well_names_to_index[name]; - WellNode& wellnode = *well_collection_.getLeafNodes()[wix]; - assert(wellnode.name() == name); - if (well_data[wix].type == PRODUCER) { - wellnode.prodSpec().guide_rate_ = lines[i].guide_rate_; - if (lines[i].phase_ == "OIL") { - wellnode.prodSpec().guide_rate_type_ = ProductionSpecification::OIL; - } else { - OPM_THROW(std::runtime_error, "Guide rate type " << lines[i].phase_ << " specified for producer " - << name << " in WGRUPCON, cannot handle."); - } - } else if (well_data[wix].type == INJECTOR) { - wellnode.injSpec().guide_rate_ = lines[i].guide_rate_; - if (lines[i].phase_ == "RAT") { - wellnode.injSpec().guide_rate_type_ = InjectionSpecification::RAT; - } else { - OPM_THROW(std::runtime_error, "Guide rate type " << lines[i].phase_ << " specified for injector " - << name << " in WGRUPCON, cannot handle."); - } - } else { - OPM_THROW(std::runtime_error, "Unknown well type " << well_data[wix].type << " for well " << name); - } - } - } - well_collection_.setWellsPointer(w_); - well_collection_.applyGroupControls(); + init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.dimensions, + grid.cell_centroids, UgGridHelpers::cell2Faces(grid), grid.face_centroids, + permeability); } - /// Destructor. WellsManager::~WellsManager() { @@ -985,141 +392,25 @@ namespace Opm well_collection_.applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase); } - void WellsManager::setupCompressedToCartesian(const UnstructuredGrid& grid, std::map& cartesian_to_compressed ) { + void WellsManager::setupCompressedToCartesian(const int* global_cell, int number_of_cells, + std::map& cartesian_to_compressed ) { // global_cell is a map from compressed cells to Cartesian grid cells. // We must make the inverse lookup. - const int* global_cell = grid.global_cell; if (global_cell) { - for (int i = 0; i < grid.number_of_cells; ++i) { + for (int i = 0; i < number_of_cells; ++i) { cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); } } else { - for (int i = 0; i < grid.number_of_cells; ++i) { + for (int i = 0; i < number_of_cells; ++i) { cartesian_to_compressed.insert(std::make_pair(i, i)); } } } - void WellsManager::createWellsFromSpecs(std::vector& wells, size_t timeStep, - const UnstructuredGrid& grid, - std::vector& well_names, - std::vector& well_data, - std::map& well_names_to_index, - const PhaseUsage& phaseUsage, - std::map cartesian_to_compressed, - const double* permeability) - { - std::vector > wellperf_data; - wellperf_data.resize(wells.size()); - int well_index = 0; - for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) { - WellConstPtr well = (*wellIter); - { // WELSPECS handling - well_names_to_index[well->name()] = well_index; - well_names.push_back(well->name()); - { - WellData wd; - // If negative (defaulted), set refdepth to a marker - // value, will be changed after getting perforation - // data to the centroid of the cell of the top well - // perforation. - wd.reference_bhp_depth = (well->getRefDepth() < 0.0) ? -1e100 : well->getRefDepth(); - wd.welspecsline = -1; - if (well->isInjector( timeStep )) - wd.type = INJECTOR; - else - wd.type = PRODUCER; - well_data.push_back(wd); - } - } - - { // COMPDAT handling - CompletionSetConstPtr completionSet = well->getCompletions(timeStep); - for (size_t c=0; csize(); c++) { - CompletionConstPtr completion = completionSet->get(c); - int i = completion->getI(); - int j = completion->getJ(); - int k = completion->getK(); - - const int* cpgdim = grid.cartdims; - int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k); - std::map::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx); - if (cgit == cartesian_to_compressed.end()) { - OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' ' - << k << " not found in grid (well = " << well->name() << ')'); - } - int cell = cgit->second; - PerfData pd; - pd.cell = cell; - if (completion->getCF() > 0.0) { - pd.well_index = completion->getCF(); - } else { - double radius = 0.5*completion->getDiameter(); - if (radius <= 0.0) { - radius = 0.5*unit::feet; - OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); - } - std::array cubical = getCubeDim(grid, cell); - const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell]; - pd.well_index = computeWellIndex(radius, cubical, cell_perm, completion->getDiameter()); - } - wellperf_data[well_index].push_back(pd); - } - } - well_index++; - } - - // Set up reference depths that were defaulted. Count perfs. - - const int num_wells = well_data.size(); - - int num_perfs = 0; - assert(grid.dimensions == 3); - for (int w = 0; w < num_wells; ++w) { - num_perfs += wellperf_data[w].size(); - if (well_data[w].reference_bhp_depth < 0.0) { - // It was defaulted. Set reference depth to minimum perforation depth. - double min_depth = 1e100; - int num_wperfs = wellperf_data[w].size(); - for (int perf = 0; perf < num_wperfs; ++perf) { - double depth = grid.cell_centroids[3*wellperf_data[w][perf].cell + 2]; - min_depth = std::min(min_depth, depth); - } - well_data[w].reference_bhp_depth = min_depth; - } - } - - // Create the well data structures. - w_ = create_wells(phaseUsage.num_phases, num_wells, num_perfs); - if (!w_) { - OPM_THROW(std::runtime_error, "Failed creating Wells struct."); - } - - - // Add wells. - for (int w = 0; w < num_wells; ++w) { - const int w_num_perf = wellperf_data[w].size(); - std::vector perf_cells(w_num_perf); - std::vector perf_prodind(w_num_perf); - for (int perf = 0; perf < w_num_perf; ++perf) { - perf_cells[perf] = wellperf_data[w][perf].cell; - perf_prodind[perf] = wellperf_data[w][perf].well_index; - } - const double* comp_frac = NULL; - // We initialize all wells with a null component fraction, - // and must (for injection wells) overwrite it later. - int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, w_num_perf, - comp_frac, &perf_cells[0], &perf_prodind[0], well_names[w].c_str(), w_); - if (!ok) { - OPM_THROW(std::runtime_error, "Failed adding well " << well_names[w] << " to Wells data structure."); - } - } - - } void WellsManager::setupWellControls(std::vector& wells, size_t timeStep, std::vector& well_names, const PhaseUsage& phaseUsage) { @@ -1137,7 +428,7 @@ namespace Opm int control_pos[5] = { -1, -1, -1, -1, -1 }; if (well->hasInjectionControl(timeStep , WellInjector::RATE)) { - control_pos[InjectionControl::RATE] = well_controls_get_num(w_->ctrls[well_index]); + control_pos[WellsManagerDetail::InjectionControl::RATE] = well_controls_get_num(w_->ctrls[well_index]); double distr[3] = { 0.0, 0.0, 0.0 }; WellInjector::TypeEnum injectorType = well->getInjectorType(timeStep); @@ -1157,7 +448,7 @@ namespace Opm } if (ok && well->hasInjectionControl(timeStep , WellInjector::RESV)) { - control_pos[InjectionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]); + control_pos[WellsManagerDetail::InjectionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]); double distr[3] = { 0.0, 0.0, 0.0 }; WellInjector::TypeEnum injectorType = well->getInjectorType(timeStep); @@ -1177,7 +468,7 @@ namespace Opm } if (ok && well->hasInjectionControl(timeStep , WellInjector::BHP)) { - control_pos[InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]); + control_pos[WellsManagerDetail::InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]); ok = append_well_controls(BHP, well->getBHPLimit(timeStep), NULL, @@ -1196,9 +487,9 @@ namespace Opm { - InjectionControl::Mode mode = InjectionControl::mode( well->getInjectorControlMode(timeStep) ); + WellsManagerDetail::InjectionControl::Mode mode = WellsManagerDetail::InjectionControl::mode( well->getInjectorControlMode(timeStep) ); int cpos = control_pos[mode]; - if (cpos == -1 && mode != InjectionControl::GRUP) { + if (cpos == -1 && mode != WellsManagerDetail::InjectionControl::GRUP) { OPM_THROW(std::runtime_error, "Control not specified in well " << well_names[well_index]); } @@ -1244,7 +535,7 @@ namespace Opm OPM_THROW(std::runtime_error, "Oil phase not active and ORAT control specified."); } - control_pos[ProductionControl::ORAT] = well_controls_get_num(w_->ctrls[well_index]); + control_pos[WellsManagerDetail::ProductionControl::ORAT] = well_controls_get_num(w_->ctrls[well_index]); double distr[3] = { 0.0, 0.0, 0.0 }; distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0; ok = append_well_controls(SURFACE_RATE, @@ -1258,7 +549,7 @@ namespace Opm if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) { OPM_THROW(std::runtime_error, "Water phase not active and WRAT control specified."); } - control_pos[ProductionControl::WRAT] = well_controls_get_num(w_->ctrls[well_index]); + control_pos[WellsManagerDetail::ProductionControl::WRAT] = well_controls_get_num(w_->ctrls[well_index]); double distr[3] = { 0.0, 0.0, 0.0 }; distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0; ok = append_well_controls(SURFACE_RATE, @@ -1272,7 +563,7 @@ namespace Opm if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) { OPM_THROW(std::runtime_error, "Gas phase not active and GRAT control specified."); } - control_pos[ProductionControl::GRAT] = well_controls_get_num(w_->ctrls[well_index]); + control_pos[WellsManagerDetail::ProductionControl::GRAT] = well_controls_get_num(w_->ctrls[well_index]); double distr[3] = { 0.0, 0.0, 0.0 }; distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0; ok = append_well_controls(SURFACE_RATE, @@ -1289,7 +580,7 @@ namespace Opm if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) { OPM_THROW(std::runtime_error, "Oil phase not active and LRAT control specified."); } - control_pos[ProductionControl::LRAT] = well_controls_get_num(w_->ctrls[well_index]); + control_pos[WellsManagerDetail::ProductionControl::LRAT] = well_controls_get_num(w_->ctrls[well_index]); double distr[3] = { 0.0, 0.0, 0.0 }; distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0; distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0; @@ -1301,7 +592,7 @@ namespace Opm } if (ok && well->hasProductionControl(timeStep , WellProducer::RESV)) { - control_pos[ProductionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]); + control_pos[WellsManagerDetail::ProductionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]); double distr[3] = { 1.0, 1.0, 1.0 }; ok = append_well_controls(RESERVOIR_RATE, -well->getResVRate(timeStep), @@ -1311,7 +602,7 @@ namespace Opm } if (ok && well->hasProductionControl(timeStep , WellProducer::BHP)) { - control_pos[ProductionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]); + control_pos[WellsManagerDetail::ProductionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]); ok = append_well_controls(BHP, well->getBHPLimit( timeStep ) , NULL, @@ -1327,9 +618,9 @@ namespace Opm OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[well_index]); } - ProductionControl::Mode mode = ProductionControl::mode(well->getProducerControlMode(timeStep)); + WellsManagerDetail::ProductionControl::Mode mode = WellsManagerDetail::ProductionControl::mode(well->getProducerControlMode(timeStep)); int cpos = control_pos[mode]; - if (cpos == -1 && mode != ProductionControl::GRUP) { + if (cpos == -1 && mode != WellsManagerDetail::ProductionControl::GRUP) { OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]); } // If it's shut, we complement the cpos diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp index 9e5c314cc..c92ba3a2c 100644 --- a/opm/core/wells/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -66,7 +66,7 @@ namespace Opm /// manage control switching does not exist. /// /// @param[in] W Existing wells object. - WellsManager(struct Wells* W); + WellsManager(struct Wells* W, bool checkCellExistence=true); /// Construct from input deck and grid. /// The permeability argument may be zero if the input contain @@ -74,13 +74,27 @@ namespace Opm /// order to approximate these by the Peaceman formula. WellsManager(const Opm::EclipseGridParser& deck, const UnstructuredGrid& grid, - const double* permeability); + const double* permeability, + bool checkCellExistence=true); + template + WellsManager(const Opm::EclipseGridParser& deck, + int num_cells, + const int* global_cell, + const int* cart_dims, + int dimensions, + CC begin_cell_centroids, + const F2C& f2c, + FC begin_face_centroids, + const double* permeability, + bool checkCellExistence=true); + WellsManager(const Opm::EclipseStateConstPtr eclipseState, const size_t timeStep, const UnstructuredGrid& grid, - const double* permeability); + const double* permeability, + bool checkCellExistence=true); /// Destructor. ~WellsManager(); @@ -134,15 +148,29 @@ namespace Opm private: + template + void init(const Opm::EclipseGridParser& deck, + int num_cells, + const int* global_cell, + const int* cart_dims, + int dimensions, + CC begin_cell_centroids, + const F2C& f2c, + FC begin_face_centroids, + const double* permeability); // Disable copying and assignment. - WellsManager(const WellsManager& other); + WellsManager(const WellsManager& other, bool checkCellExistence=true); WellsManager& operator=(const WellsManager& other); - static void setupCompressedToCartesian(const UnstructuredGrid& grid, std::map& cartesian_to_compressed ); + static void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map& cartesian_to_compressed ); void setupWellControls(std::vector& wells, size_t timeStep, std::vector& well_names, const PhaseUsage& phaseUsage); - + template void createWellsFromSpecs( std::vector& wells, size_t timeStep, - const UnstructuredGrid& grid, + const C2F& c2f, + const int* cart_dims, + FC begin_face_centroids, + CC begin_cell_centroids, + int dimensions, std::vector& well_names, std::vector& well_data, std::map & well_names_to_index, @@ -157,9 +185,10 @@ namespace Opm // Data Wells* w_; WellCollection well_collection_; + bool checkCellExistence_; }; } // namespace Opm - +#include "WellsManager_impl.hpp" #endif // OPM_WELLSMANAGER_HEADER_INCLUDED diff --git a/opm/core/wells/WellsManager_impl.hpp b/opm/core/wells/WellsManager_impl.hpp new file mode 100644 index 000000000..203b10c25 --- /dev/null +++ b/opm/core/wells/WellsManager_impl.hpp @@ -0,0 +1,819 @@ +#include +#include + +#include + +namespace WellsManagerDetail +{ + + + namespace ProductionControl + { + enum Mode { ORAT, WRAT, GRAT, + LRAT, CRAT, RESV, + BHP , THP , GRUP }; + /* + namespace Details { + std::map + init_mode_map(); + } // namespace Details + */ + Mode mode(const std::string& control); + + + Mode mode(Opm::WellProducer::ControlModeEnum controlMode); + } // namespace ProductionControl + + + namespace InjectionControl + { + enum Mode { RATE, RESV, BHP, + THP, GRUP }; + /* + namespace Details { + std::map + init_mode_map(); + } // namespace Details + */ + Mode mode(const std::string& control); + + Mode mode(Opm::WellInjector::ControlModeEnum controlMode); + + } // namespace InjectionControl + +double computeWellIndex(const double radius, + const std::array& cubical, + const double* cell_permeability, + const double skin_factor); + +template +std::array getCubeDim(const C2F& c2f, FC begin_face_centroids, int dimensions, + int cell) +{ + using namespace std; + std::array cube; + //typedef Opm::UgGridHelpers::Cell2FacesTraits::Type Cell2Faces; + //Cell2Faces c2f=Opm::UgGridHelpers::cell2Faces(grid); + typedef typename C2F::row_type FaceRow; + FaceRow faces=c2f[cell]; + typename FaceRow::const_iterator face=faces.begin(); + int num_local_faces = faces.end()-face; + vector x(num_local_faces); + vector y(num_local_faces); + vector z(num_local_faces); + for (int lf=0; lf +void WellsManager::createWellsFromSpecs(std::vector& wells, size_t timeStep, + const C2F& c2f, + const int* cart_dims, + FC begin_face_centroids, + CC begin_cell_centroids, + int dimensions, + std::vector& well_names, + std::vector& well_data, + std::map& well_names_to_index, + const PhaseUsage& phaseUsage, + std::map cartesian_to_compressed, + const double* permeability) +{ + std::vector > wellperf_data; + wellperf_data.resize(wells.size()); + + int well_index = 0; + for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) { + WellConstPtr well = (*wellIter); + { // WELSPECS handling + well_names_to_index[well->name()] = well_index; + well_names.push_back(well->name()); + { + WellData wd; + // If negative (defaulted), set refdepth to a marker + // value, will be changed after getting perforation + // data to the centroid of the cell of the top well + // perforation. + wd.reference_bhp_depth = (well->getRefDepth() < 0.0) ? -1e100 : well->getRefDepth(); + wd.welspecsline = -1; + if (well->isInjector( timeStep )) + wd.type = INJECTOR; + else + wd.type = PRODUCER; + well_data.push_back(wd); + } + } + + { // COMPDAT handling + CompletionSetConstPtr completionSet = well->getCompletions(timeStep); + for (size_t c=0; csize(); c++) { + CompletionConstPtr completion = completionSet->get(c); + int i = completion->getI(); + int j = completion->getJ(); + int k = completion->getK(); + + const int* cpgdim = cart_dims; + int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k); + std::map::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx); + if (cgit == cartesian_to_compressed.end()) { + if (checkCellExistence_) + OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' ' + << k << " not found in grid (well = " << well->name() << ')'); + continue; + } + int cell = cgit->second; + PerfData pd; + pd.cell = cell; + if (completion->getCF() > 0.0) { + pd.well_index = completion->getCF(); + } else { + double radius = 0.5*completion->getDiameter(); + if (radius <= 0.0) { + radius = 0.5*unit::feet; + OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); + } + std::array cubical = WellsManagerDetail::getCubeDim(c2f, begin_face_centroids, + dimensions, cell); + const double* cell_perm = &permeability[dimensions*dimensions*cell]; + pd.well_index = WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm, completion->getDiameter()); + } + wellperf_data[well_index].push_back(pd); + } + } + well_index++; + } + + // Set up reference depths that were defaulted. Count perfs. + + const int num_wells = well_data.size(); + + int num_perfs = 0; + assert(dimensions == 3); + for (int w = 0; w < num_wells; ++w) { + num_perfs += wellperf_data[w].size(); + if (well_data[w].reference_bhp_depth < 0.0) { + // It was defaulted. Set reference depth to minimum perforation depth. + double min_depth = 1e100; + int num_wperfs = wellperf_data[w].size(); + for (int perf = 0; perf < num_wperfs; ++perf) { + double depth = UgGridHelpers + ::getCoordinate(UgGridHelpers::increment(begin_cell_centroids, + wellperf_data[w][perf].cell, + dimensions), + 2); + min_depth = std::min(min_depth, depth); + } + well_data[w].reference_bhp_depth = min_depth; + } + } + + // Create the well data structures. + w_ = create_wells(phaseUsage.num_phases, num_wells, num_perfs); + if (!w_) { + OPM_THROW(std::runtime_error, "Failed creating Wells struct."); + } + + + // Add wells. + for (int w = 0; w < num_wells; ++w) { + const int w_num_perf = wellperf_data[w].size(); + std::vector perf_cells(w_num_perf); + std::vector perf_prodind(w_num_perf); + for (int perf = 0; perf < w_num_perf; ++perf) { + perf_cells[perf] = wellperf_data[w][perf].cell; + perf_prodind[perf] = wellperf_data[w][perf].well_index; + } + const double* comp_frac = NULL; + // We initialize all wells with a null component fraction, + // and must (for injection wells) overwrite it later. + int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, w_num_perf, + comp_frac, &perf_cells[0], &perf_prodind[0], well_names[w].c_str(), w_); + if (!ok) { + OPM_THROW(std::runtime_error, "Failed adding well " << well_names[w] << " to Wells data structure."); + } + } +} + +template +WellsManager::WellsManager(const Opm::EclipseGridParser& deck, + int number_of_cells, + const int* global_cell, + const int* cart_dims, + int dimensions, + CC begin_cell_centroids, + const C2F& c2f, + FC begin_face_centroids, + const double* permeability, + bool checkCellExistence) + : w_(0), checkCellExistence_(checkCellExistence) +{ + init(deck, number_of_cells, global_cell, cart_dims, dimensions, + begin_cell_centroids, c2f, begin_face_centroids, permeability); +} + + /// Construct wells from deck. +template +void WellsManager::init(const Opm::EclipseGridParser& deck, + int number_of_cells, + const int* global_cell, + const int* cart_dims, + int dimensions, + CC begin_cell_centroids, + const C2F& c2f, + FC begin_face_centroids, + const double* permeability) +{ + if (dimensions != 3) { + OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional."); + } + // NOTE: Implementation copied and modified from dune-porsol's class BlackoilWells. + std::vector keywords; + keywords.push_back("WELSPECS"); + keywords.push_back("COMPDAT"); + // keywords.push_back("WELTARG"); + if (!deck.hasFields(keywords)) { + OPM_MESSAGE("Missing well keywords in deck, initializing no wells."); + return; + } + if (!(deck.hasField("WCONINJE") || deck.hasField("WCONPROD")) ) { + OPM_THROW(std::runtime_error, "Needed field is missing in file"); + } + + // Obtain phase usage data. + PhaseUsage pu = phaseUsageFromDeck(deck); + + // These data structures will be filled in this constructor, + // then used to initialize the Wells struct. + std::vector well_names; + std::vector well_data; + std::vector > wellperf_data; + + // For easy lookup: + std::map well_names_to_index; + typedef std::map::const_iterator WNameIt; + + // Get WELSPECS data. + // It is allowed to have multiple lines corresponding to + // the same well, in which case the last one encountered + // is the one used. + const WELSPECS& welspecs = deck.getWELSPECS(); + const int num_welspecs = welspecs.welspecs.size(); + well_names.reserve(num_welspecs); + well_data.reserve(num_welspecs); + for (int w = 0; w < num_welspecs; ++w) { + // First check if this well has already been encountered. + // If so, we modify it's data instead of appending a new well + // to the well_data and well_names vectors. + const std::string& name = welspecs.welspecs[w].name_; + const double refdepth = welspecs.welspecs[w].datum_depth_BHP_; + WNameIt wit = well_names_to_index.find(name); + if (wit == well_names_to_index.end()) { + // New well, append data. + well_names_to_index[welspecs.welspecs[w].name_] = well_data.size(); + well_names.push_back(name); + WellData wd; + // If negative (defaulted), set refdepth to a marker + // value, will be changed after getting perforation + // data to the centroid of the cell of the top well + // perforation. + wd.reference_bhp_depth = (refdepth < 0.0) ? -1e100 : refdepth; + wd.welspecsline = w; + well_data.push_back(wd); + } else { + // Existing well, change data. + const int wix = wit->second; + well_data[wix].reference_bhp_depth = (refdepth < 0.0) ? -1e100 : refdepth; + well_data[wix].welspecsline = w; + } + } + const int num_wells = well_data.size(); + wellperf_data.resize(num_wells); + + + // global_cell is a map from compressed cells to Cartesian grid cells. + // We must make the inverse lookup. + const int* cpgdim = cart_dims; + std::map cartesian_to_compressed; + + if (global_cell) { + for (int i = 0; i < number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); + } + } + else { + for (int i = 0; i < number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(i, i)); + } + } + + // Get COMPDAT data + // It is *not* allowed to have multiple lines corresponding to + // the same perforation! + const COMPDAT& compdat = deck.getCOMPDAT(); + const int num_compdat = compdat.compdat.size(); + for (int kw = 0; kw < num_compdat; ++kw) { + // Extract well name, or the part of the well name that + // comes before the '*'. + std::string name = compdat.compdat[kw].well_; + std::string::size_type len = name.find('*'); + if (len != std::string::npos) { + name = name.substr(0, len); + } + // Look for well with matching name. + bool found = false; + for (int wix = 0; wix < num_wells; ++wix) { + if (well_names[wix].compare(0,len, name) == 0) { // equal + // Extract corresponding WELSPECS defintion for + // purpose of default location specification. + const WelspecsLine& wspec = welspecs.welspecs[well_data[wix].welspecsline]; + + // We have a matching name. + int ix = compdat.compdat[kw].grid_ind_[0] - 1; + int jy = compdat.compdat[kw].grid_ind_[1] - 1; + int kz1 = compdat.compdat[kw].grid_ind_[2] - 1; + int kz2 = compdat.compdat[kw].grid_ind_[3] - 1; + + if (ix < 0) { + // Defaulted I location. Extract from WELSPECS. + ix = wspec.I_ - 1; + } + if (jy < 0) { + // Defaulted J location. Extract from WELSPECS. + jy = wspec.J_ - 1; + } + if (kz1 < 0) { + // Defaulted KZ1. Use top layer. + kz1 = 0; + } + if (kz2 < 0) { + // Defaulted KZ2. Use bottom layer. + kz2 = cpgdim[2] - 1; + } + + for (int kz = kz1; kz <= kz2; ++kz) { + int cart_grid_indx = ix + cpgdim[0]*(jy + cpgdim[1]*kz); + std::map::const_iterator cgit = + cartesian_to_compressed.find(cart_grid_indx); + if (cgit == cartesian_to_compressed.end()) { + if(checkCellExistence_) + OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << ix << ' ' << jy << ' ' + << kz << " not found in grid (well = " << name << ')'); + continue; + } + int cell = cgit->second; + PerfData pd; + pd.cell = cell; + if (compdat.compdat[kw].connect_trans_fac_ > 0.0) { + pd.well_index = compdat.compdat[kw].connect_trans_fac_; + } else { + double radius = 0.5*compdat.compdat[kw].diameter_; + if (radius <= 0.0) { + radius = 0.5*unit::feet; + OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); + } + std::array cubical = WellsManagerDetail::getCubeDim(c2f, + begin_face_centroids, + dimensions, + cell); + const double* cell_perm = &permeability[dimensions*dimensions*cell]; + pd.well_index = WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm, + compdat.compdat[kw].skin_factor_); + } + wellperf_data[wix].push_back(pd); + } + found = true; + break; + } + } + if (!found) { + OPM_THROW(std::runtime_error, "Undefined well name: " << compdat.compdat[kw].well_ + << " in COMPDAT"); + } + } + + // Set up reference depths that were defaulted. Count perfs. + int num_perfs = 0; + assert(dimensions == 3); + for (int w = 0; w < num_wells; ++w) { + num_perfs += wellperf_data[w].size(); + if (well_data[w].reference_bhp_depth < 0.0) { + // It was defaulted. Set reference depth to minimum perforation depth. + double min_depth = 1e100; + int num_wperfs = wellperf_data[w].size(); + for (int perf = 0; perf < num_wperfs; ++perf) { + double depth = UgGridHelpers:: + getCoordinate(UgGridHelpers::increment(begin_cell_centroids, + wellperf_data[w][perf].cell, + dimensions), 2); + min_depth = std::min(min_depth, depth); + } + well_data[w].reference_bhp_depth = min_depth; + } + } + + // Create the well data structures. + w_ = create_wells(pu.num_phases, num_wells, num_perfs); + if (!w_) { + OPM_THROW(std::runtime_error, "Failed creating Wells struct."); + } + + // Classify wells + if (deck.hasField("WCONINJE")) { + const std::vector& lines = deck.getWCONINJE().wconinje; + for (size_t i = 0 ; i < lines.size(); ++i) { + const std::map::const_iterator it = well_names_to_index.find(lines[i].well_); + if (it != well_names_to_index.end()) { + const int well_index = it->second; + well_data[well_index].type = INJECTOR; + } else { + OPM_THROW(std::runtime_error, "Unseen well name: " << lines[i].well_ << " first seen in WCONINJE"); + } + } + } + if (deck.hasField("WCONPROD")) { + const std::vector& lines = deck.getWCONPROD().wconprod; + for (size_t i = 0; i < lines.size(); ++i) { + const std::map::const_iterator it = well_names_to_index.find(lines[i].well_); + if (it != well_names_to_index.end()) { + const int well_index = it->second; + well_data[well_index].type = PRODUCER; + } else { + OPM_THROW(std::runtime_error, "Unseen well name: " << lines[i].well_ << " first seen in WCONPROD"); + } + + } + } + + // Add wells. + for (int w = 0; w < num_wells; ++w) { + const int w_num_perf = wellperf_data[w].size(); + std::vector perf_cells(w_num_perf); + std::vector perf_prodind(w_num_perf); + for (int perf = 0; perf < w_num_perf; ++perf) { + perf_cells[perf] = wellperf_data[w][perf].cell; + perf_prodind[perf] = wellperf_data[w][perf].well_index; + } + const double* comp_frac = NULL; + // We initialize all wells with a null component fraction, + // and must (for injection wells) overwrite it later. + int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, w_num_perf, + comp_frac, &perf_cells[0], &perf_prodind[0], well_names[w].c_str(), w_); + if (!ok) { + OPM_THROW(std::runtime_error, "Failed adding well " << well_names[w] << " to Wells data structure."); + } + } + + // Get WCONINJE data, add injection controls to wells. + // It is allowed to have multiple lines corresponding to + // the same well, in which case the last one encountered + // is the one used. + if (deck.hasField("WCONINJE")) { + const WCONINJE& wconinjes = deck.getWCONINJE(); + const int num_wconinjes = wconinjes.wconinje.size(); + for (int kw = 0; kw < num_wconinjes; ++kw) { + const WconinjeLine& wci_line = wconinjes.wconinje[kw]; + // Extract well name, or the part of the well name that + // comes before the '*'. + std::string name = wci_line.well_; + std::string::size_type len = name.find('*'); + if (len != std::string::npos) { + name = name.substr(0, len); + } + bool well_found = false; + for (int wix = 0; wix < num_wells; ++wix) { + if (well_names[wix].compare(0,len, name) == 0) { //equal + well_found = true; + assert(well_data[wix].type == w_->type[wix]); + if (well_data[wix].type != INJECTOR) { + OPM_THROW(std::runtime_error, "Found WCONINJE entry for a non-injector well: " << well_names[wix]); + } + + // Add all controls that are present in well. + // First we must clear existing controls, in case the + // current WCONINJE line is modifying earlier controls. + clear_well_controls(wix, w_); + int ok = 1; + int control_pos[5] = { -1, -1, -1, -1, -1 }; + if (ok && wci_line.surface_flow_max_rate_ >= 0.0) { + control_pos[WellsManagerDetail::InjectionControl::RATE] = well_controls_get_num(w_->ctrls[wix]); + double distr[3] = { 0.0, 0.0, 0.0 }; + if (wci_line.injector_type_ == "WATER") { + distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; + } else if (wci_line.injector_type_ == "OIL") { + distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; + } else if (wci_line.injector_type_ == "GAS") { + distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; + } else { + OPM_THROW(std::runtime_error, "Injector type " << wci_line.injector_type_ << " not supported." + "WellsManager only supports WATER, OIL and GAS injector types."); + } + ok = append_well_controls(SURFACE_RATE, wci_line.surface_flow_max_rate_, + distr, wix, w_); + } + if (ok && wci_line.reservoir_flow_max_rate_ >= 0.0) { + control_pos[WellsManagerDetail::InjectionControl::RESV] = well_controls_get_num(w_->ctrls[wix]); + double distr[3] = { 0.0, 0.0, 0.0 }; + if (wci_line.injector_type_ == "WATER") { + distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; + } else if (wci_line.injector_type_ == "OIL") { + distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; + } else if (wci_line.injector_type_ == "GAS") { + distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; + } else { + OPM_THROW(std::runtime_error, "Injector type " << wci_line.injector_type_ << " not supported." + "WellsManager only supports WATER, OIL and GAS injector types."); + } + ok = append_well_controls(RESERVOIR_RATE, wci_line.reservoir_flow_max_rate_, + distr, wix, w_); + } + if (ok && wci_line.BHP_limit_ > 0.0) { + control_pos[WellsManagerDetail::InjectionControl::BHP] = well_controls_get_num(w_->ctrls[wix]); + ok = append_well_controls(BHP, wci_line.BHP_limit_, + NULL, wix, w_); + } + if (ok && wci_line.THP_limit_ > 0.0) { + OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[wix]); + } + if (!ok) { + OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[wix]); + } + WellsManagerDetail::InjectionControl::Mode mode = WellsManagerDetail::InjectionControl::mode(wci_line.control_mode_); + int cpos = control_pos[mode]; + if (cpos == -1 && mode != WellsManagerDetail::InjectionControl::GRUP) { + OPM_THROW(std::runtime_error, "Control for " << wci_line.control_mode_ << " not specified in well " << well_names[wix]); + } + // We need to check if the well is shut or not + if (wci_line.open_shut_flag_ == "SHUT") { + cpos = ~cpos; + } + set_current_control(wix, cpos, w_); + + // Set well component fraction. + double cf[3] = { 0.0, 0.0, 0.0 }; + if (wci_line.injector_type_[0] == 'W') { + if (!pu.phase_used[BlackoilPhases::Aqua]) { + OPM_THROW(std::runtime_error, "Water phase not used, yet found water-injecting well."); + } + cf[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; + } else if (wci_line.injector_type_[0] == 'O') { + if (!pu.phase_used[BlackoilPhases::Liquid]) { + OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-injecting well."); + } + cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; + } else if (wci_line.injector_type_[0] == 'G') { + if (!pu.phase_used[BlackoilPhases::Vapour]) { + OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-injecting well."); + } + cf[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; + } + std::copy(cf, cf + pu.num_phases, w_->comp_frac + wix*pu.num_phases); + } + } + if (!well_found) { + OPM_THROW(std::runtime_error, "Undefined well name: " << wci_line.well_ + << " in WCONINJE"); + } + } + } + + // Get WCONPROD data + // It is allowed to have multiple lines corresponding to + // the same well, in which case the last one encountered + // is the one used. + if (deck.hasField("WCONPROD")) { + const WCONPROD& wconprods = deck.getWCONPROD(); + const int num_wconprods = wconprods.wconprod.size(); + for (int kw = 0; kw < num_wconprods; ++kw) { + const WconprodLine& wcp_line = wconprods.wconprod[kw]; + std::string name = wcp_line.well_; + std::string::size_type len = name.find('*'); + if (len != std::string::npos) { + name = name.substr(0, len); + } + bool well_found = false; + for (int wix = 0; wix < num_wells; ++wix) { + if (well_names[wix].compare(0,len, name) == 0) { //equal + well_found = true; + assert(well_data[wix].type == w_->type[wix]); + if (well_data[wix].type != PRODUCER) { + OPM_THROW(std::runtime_error, "Found WCONPROD entry for a non-producer well: " << well_names[wix]); + } + // Add all controls that are present in well. + // First we must clear existing controls, in case the + // current WCONPROD line is modifying earlier controls. + clear_well_controls(wix, w_); + int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 }; + int ok = 1; + if (ok && wcp_line.oil_max_rate_ >= 0.0) { + if (!pu.phase_used[BlackoilPhases::Liquid]) { + OPM_THROW(std::runtime_error, "Oil phase not active and ORAT control specified."); + } + control_pos[WellsManagerDetail::ProductionControl::ORAT] = well_controls_get_num(w_->ctrls[wix]); + double distr[3] = { 0.0, 0.0, 0.0 }; + distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; + ok = append_well_controls(SURFACE_RATE, -wcp_line.oil_max_rate_, + distr, wix, w_); + } + if (ok && wcp_line.water_max_rate_ >= 0.0) { + if (!pu.phase_used[BlackoilPhases::Aqua]) { + OPM_THROW(std::runtime_error, "Water phase not active and WRAT control specified."); + } + control_pos[WellsManagerDetail::ProductionControl::WRAT] = well_controls_get_num(w_->ctrls[wix]); + double distr[3] = { 0.0, 0.0, 0.0 }; + distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; + ok = append_well_controls(SURFACE_RATE, -wcp_line.water_max_rate_, + distr, wix, w_); + } + if (ok && wcp_line.gas_max_rate_ >= 0.0) { + if (!pu.phase_used[BlackoilPhases::Vapour]) { + OPM_THROW(std::runtime_error, "Gas phase not active and GRAT control specified."); + } + control_pos[WellsManagerDetail::ProductionControl::GRAT] = well_controls_get_num(w_->ctrls[wix]); + double distr[3] = { 0.0, 0.0, 0.0 }; + distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; + ok = append_well_controls(SURFACE_RATE, -wcp_line.gas_max_rate_, + distr, wix, w_); + } + if (ok && wcp_line.liquid_max_rate_ >= 0.0) { + if (!pu.phase_used[BlackoilPhases::Aqua]) { + OPM_THROW(std::runtime_error, "Water phase not active and LRAT control specified."); + } + if (!pu.phase_used[BlackoilPhases::Liquid]) { + OPM_THROW(std::runtime_error, "Oil phase not active and LRAT control specified."); + } + control_pos[WellsManagerDetail::ProductionControl::LRAT] = well_controls_get_num(w_->ctrls[wix]); + double distr[3] = { 0.0, 0.0, 0.0 }; + distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; + distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; + ok = append_well_controls(SURFACE_RATE, -wcp_line.liquid_max_rate_, + distr, wix, w_); + } + if (ok && wcp_line.reservoir_flow_max_rate_ >= 0.0) { + control_pos[WellsManagerDetail::ProductionControl::RESV] = well_controls_get_num(w_->ctrls[wix]); + double distr[3] = { 1.0, 1.0, 1.0 }; + ok = append_well_controls(RESERVOIR_RATE, -wcp_line.reservoir_flow_max_rate_, + distr, wix, w_); + } + if (ok && wcp_line.BHP_limit_ > 0.0) { + control_pos[WellsManagerDetail::ProductionControl::BHP] = well_controls_get_num(w_->ctrls[wix]); + ok = append_well_controls(BHP, wcp_line.BHP_limit_, + NULL, wix, w_); + } + if (ok && wcp_line.THP_limit_ > 0.0) { + OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[wix]); + } + if (!ok) { + OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[wix]); + } + WellsManagerDetail::ProductionControl::Mode mode = WellsManagerDetail::ProductionControl::mode(wcp_line.control_mode_); + int cpos = control_pos[mode]; + if (cpos == -1 && mode != WellsManagerDetail::ProductionControl::GRUP) { + OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[wix]); + } + // If it's shut, we complement the cpos + if (wcp_line.open_shut_flag_ == "SHUT") { + cpos = ~cpos; // So we can easily retrieve the cpos later + } + set_current_control(wix, cpos, w_); + } + } + if (!well_found) { + OPM_THROW(std::runtime_error, "Undefined well name: " << wcp_line.well_ + << " in WCONPROD"); + } + } + } + + // Get WELTARG data + if (deck.hasField("WELTARG")) { + OPM_THROW(std::runtime_error, "We currently do not handle WELTARG."); + /* + const WELTARG& weltargs = deck.getWELTARG(); + const int num_weltargs = weltargs.weltarg.size(); + for (int kw = 0; kw < num_weltargs; ++kw) { + std::string name = weltargs.weltarg[kw].well_; + std::string::size_type len = name.find('*'); + if (len != std::string::npos) { + name = name.substr(0, len); + } + bool well_found = false; + for (int wix = 0; wix < num_wells; ++wix) { + if (well_names[wix].compare(0,len, name) == 0) { //equal + well_found = true; + well_data[wix].target = weltargs.weltarg[kw].new_value_; + break; + } + } + if (!well_found) { + OPM_THROW(std::runtime_error, "Undefined well name: " << weltargs.weltarg[kw].well_ + << " in WELTARG"); + } + } + */ + } + + // Debug output. +#define EXTRA_OUTPUT +#ifdef EXTRA_OUTPUT + /* + std::cout << "\t WELL DATA" << std::endl; + for(int i = 0; i< num_wells; ++i) { + std::cout << i << ": " << well_data[i].type << " " + << well_data[i].control << " " << well_data[i].target + << std::endl; + } + + std::cout << "\n\t PERF DATA" << std::endl; + for(int i=0; i< int(wellperf_data.size()); ++i) { + for(int j=0; j< int(wellperf_data[i].size()); ++j) { + std::cout << i << ": " << wellperf_data[i][j].cell << " " + << wellperf_data[i][j].well_index << std::endl; + } + } + */ +#endif + + if (deck.hasField("WELOPEN")) { + const WELOPEN& welopen = deck.getWELOPEN(); + for (size_t i = 0; i < welopen.welopen.size(); ++i) { + WelopenLine line = welopen.welopen[i]; + std::string wellname = line.well_; + std::map::const_iterator it = well_names_to_index.find(wellname); + if (it == well_names_to_index.end()) { + OPM_THROW(std::runtime_error, "Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS."); + } + const int index = it->second; + if (line.openshutflag_ == "SHUT") { + well_controls_shut_well( w_->ctrls[index] ); + } else if (line.openshutflag_ == "OPEN") { + well_controls_open_well( w_->ctrls[index] ); + } else { + OPM_THROW(std::runtime_error, "Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT."); + } + } + } + + // Build the well_collection_ well group hierarchy. + if (deck.hasField("GRUPTREE")) { + std::cout << "Found gruptree" << std::endl; + const GRUPTREE& gruptree = deck.getGRUPTREE(); + std::map::const_iterator it = gruptree.tree.begin(); + for( ; it != gruptree.tree.end(); ++it) { + well_collection_.addChild(it->first, it->second, deck); + } + } + for (size_t i = 0; i < welspecs.welspecs.size(); ++i) { + WelspecsLine line = welspecs.welspecs[i]; + well_collection_.addChild(line.name_, line.group_, deck); + } + + + + // Set the guide rates: + if (deck.hasField("WGRUPCON")) { + std::cout << "Found Wgrupcon" << std::endl; + WGRUPCON wgrupcon = deck.getWGRUPCON(); + const std::vector& lines = wgrupcon.wgrupcon; + std::cout << well_collection_.getLeafNodes().size() << std::endl; + for (size_t i = 0; i < lines.size(); i++) { + std::string name = lines[i].well_; + const int wix = well_names_to_index[name]; + WellNode& wellnode = *well_collection_.getLeafNodes()[wix]; + assert(wellnode.name() == name); + if (well_data[wix].type == PRODUCER) { + wellnode.prodSpec().guide_rate_ = lines[i].guide_rate_; + if (lines[i].phase_ == "OIL") { + wellnode.prodSpec().guide_rate_type_ = ProductionSpecification::OIL; + } else { + OPM_THROW(std::runtime_error, "Guide rate type " << lines[i].phase_ << " specified for producer " + << name << " in WGRUPCON, cannot handle."); + } + } else if (well_data[wix].type == INJECTOR) { + wellnode.injSpec().guide_rate_ = lines[i].guide_rate_; + if (lines[i].phase_ == "RAT") { + wellnode.injSpec().guide_rate_type_ = InjectionSpecification::RAT; + } else { + OPM_THROW(std::runtime_error, "Guide rate type " << lines[i].phase_ << " specified for injector " + << name << " in WGRUPCON, cannot handle."); + } + } else { + OPM_THROW(std::runtime_error, "Unknown well type " << well_data[wix].type << " for well " << name); + } + } + } + well_collection_.setWellsPointer(w_); + well_collection_.applyGroupControls(); +} + +} // end namespace Opm