From d5f470cb686d08b4080022ae8b12be61708e154d Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Mon, 17 Feb 2014 13:23:01 +0100 Subject: [PATCH 1/8] Refactored parts needed for Blackoil in autodiff to get rid of UG dependency. This patch refactors (hopefully) all parts of opm-core that are needed by the fully implicite black oil solver in opm-autodiff and that inherently relied on UnstructuredGrid. We added a new simple grid interface consisting out of free functions that will allow us to use CpGrid without copying it to an UnstructuredGrid by the means of the GridAdapter. Using this interface we have add methods that allow specifying the grid information (global_cell, cartdims, etc.) wherever possible to prevent introducing grid parameters for the type of the grid. Unfortunately this was not possible everywhere. --- opm/core/props/BlackoilPropertiesFromDeck.cpp | 75 +--- opm/core/props/BlackoilPropertiesFromDeck.hpp | 19 + opm/core/props/rock/RockFromDeck.cpp | 40 ++- opm/core/props/rock/RockFromDeck.hpp | 17 +- .../props/satfunc/SaturationPropsFromDeck.hpp | 33 +- .../satfunc/SaturationPropsFromDeck_impl.hpp | 124 +++++-- opm/core/simulator/BlackoilState.cpp | 13 +- opm/core/simulator/BlackoilState.hpp | 4 +- opm/core/simulator/initState.hpp | 111 +++++- opm/core/simulator/initState_impl.hpp | 338 ++++++++++++------ opm/core/utility/miscUtilities.cpp | 28 +- opm/core/utility/miscUtilities.hpp | 21 ++ opm/core/utility/miscUtilities_impl.hpp | 41 +++ opm/core/wells/WellsManager.cpp | 78 ++-- opm/core/wells/WellsManager.hpp | 13 +- 15 files changed, 667 insertions(+), 288 deletions(-) create mode 100644 opm/core/utility/miscUtilities_impl.hpp diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp index e476f015f..e5838607a 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp @@ -28,83 +28,18 @@ namespace Opm const UnstructuredGrid& grid, bool init_rock) { - if (init_rock){ - rock_.init(deck, grid); - } - pvt_.init(deck, 200); - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); - satprops_.reset(ptr); - ptr->init(deck, grid, 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(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids, + grid.dimensions, init_rock); } + BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, const UnstructuredGrid& grid, const parameter::ParameterGroup& param, bool init_rock) { - if(init_rock){ - rock_.init(deck, grid); - } - - const int pvt_samples = param.getDefault("pvt_tab_size", 200); - pvt_.init(deck, 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 (deck.hasField("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(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() diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index 49072b29e..803ce3bd0 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -184,6 +184,24 @@ 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); + RockFromDeck rock_; BlackoilPvtProperties pvt_; std::unique_ptr satprops_; @@ -197,5 +215,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 730363434..3bbe2b9fa 100644 --- a/opm/core/props/rock/RockFromDeck.cpp +++ b/opm/core/props/rock/RockFromDeck.cpp @@ -58,22 +58,34 @@ 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::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]; } } @@ -82,16 +94,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); @@ -113,13 +126,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) { diff --git a/opm/core/props/rock/RockFromDeck.hpp b/opm/core/props/rock/RockFromDeck.hpp index e277ee2a5..8e6ee01e1 100644 --- a/opm/core/props/rock/RockFromDeck.hpp +++ b/opm/core/props/rock/RockFromDeck.hpp @@ -43,6 +43,16 @@ namespace Opm void init(const EclipseGridParser& deck, const UnstructuredGrid& grid); + /// 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); + /// \return D, the number of spatial dimensions. Always 3 for deck input. int numDimensions() const { @@ -71,9 +81,12 @@ namespace Opm private: void assignPorosity(const EclipseGridParser& parser, - 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); std::vector porosity_; diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp index f2892e12a..a9177fab1 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp @@ -60,6 +60,28 @@ namespace Opm const UnstructuredGrid& grid, const int samples); + /// Initialize from deck and grid. + /// \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); + /// \return P, the number of phases. int numPhases() const; @@ -123,11 +145,14 @@ namespace Opm typedef SatFuncSet Funcs; const Funcs& funcForCell(const int cell) const; - + template void initEPS(const EclipseGridParser& deck, - const UnstructuredGrid& grid, - const std::string& keyword, - std::vector& scaleparam); + int number_of_cells, + const int* global_cell, + const T& begin_cell_centroids, + int dimensions, + const std::string& keyword, + std::vector& scaleparam); void relpermEPS(const double *s, const int cell, double *kr, double *dkrds= 0) const; }; diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp index 0b5d61791..5951aaa50 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp @@ -46,6 +46,19 @@ 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); @@ -61,11 +74,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; } } @@ -118,14 +129,22 @@ namespace Opm } } do_eps_ = true; - initEPS(deck, grid, std::string("SWCR"), eps_.swcr_); - initEPS(deck, grid, std::string("SWL"), eps_.swl_); - initEPS(deck, grid, std::string("SWU"), eps_.swu_); - initEPS(deck, grid, std::string("SOWCR"), eps_.sowcr_); - initEPS(deck, grid, std::string("KRW"), eps_.krw_); - initEPS(deck, grid, std::string("KRWR"), eps_.krwr_); - initEPS(deck, grid, std::string("KRO"), eps_.kro_); - initEPS(deck, grid, std::string("KRORW"), eps_.krorw_); + initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions, + std::string("SWCR"), eps_.swcr_); + initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions, + std::string("SWL"), eps_.swl_); + initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions, + std::string("SWU"), eps_.swu_); + initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions, + std::string("SOWCR"), eps_.sowcr_); + initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions, + std::string("KRW"), eps_.krw_); + initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions, + std::string("KRWR"), eps_.krwr_); + initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions, + std::string("KRO"), eps_.kro_); + initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions, + std::string("KRORW"), eps_.krorw_); } } @@ -253,10 +272,45 @@ namespace Opm return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]]; } +namespace +{ + + template + const T* increment(T* cc, int i, int dim) + { + return cc+(i*dim); + } + + template + T increment(const T& t, int i, int) + { + return t+i; + } + + template + double getCoordinate(T* cc, int i) + { + return cc[i]; + } + + template + double getCoordinate(T t, int i) + { + return t->center()[i]; + } + + +} + + // 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, const std::string& keyword, std::vector& scaleparam) { @@ -273,29 +327,29 @@ namespace Opm if (keyword == std::string("SWL")) { if (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 { @@ -358,14 +411,13 @@ 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 = getCoordinate(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); } 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..5841592a6 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,27 @@ 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, + 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); } // namespace Opm #include diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index 12dbe7b52..bfafd5d6b 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 @@ -47,17 +48,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 { @@ -74,8 +79,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, @@ -91,10 +98,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); @@ -113,8 +122,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, @@ -123,14 +134,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; } @@ -139,8 +150,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, @@ -149,7 +162,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); } @@ -181,8 +195,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, @@ -193,11 +209,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); } @@ -266,30 +282,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]); @@ -318,11 +346,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); @@ -336,11 +385,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); } @@ -353,30 +400,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"); @@ -388,20 +439,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 = 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 = 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); } @@ -412,13 +466,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); @@ -431,11 +505,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); } @@ -448,26 +520,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); } @@ -478,6 +554,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); @@ -485,7 +581,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."); @@ -499,17 +595,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 @@ -519,8 +616,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]; @@ -533,8 +630,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]; @@ -550,8 +647,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]; @@ -565,7 +662,8 @@ namespace Opm } // Finally, init face pressures. - initFacePressure(grid, state); + initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids, + begin_cell_centroids, state); } @@ -579,10 +677,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) { @@ -616,6 +722,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."); } @@ -627,27 +742,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]; } @@ -657,7 +771,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; @@ -667,11 +781,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) @@ -689,10 +803,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) @@ -710,9 +824,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]; @@ -744,24 +858,44 @@ namespace Opm const double gravity, State& state) { - initStateFromDeck(grid, props, deck, gravity, state); + initBlackoilStateFromDeck(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, deck, gravity, state); + } + + template + void initBlackoilStateFromDeck(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) + { + 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..73dd738ff 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -401,24 +401,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 diff --git a/opm/core/utility/miscUtilities.hpp b/opm/core/utility/miscUtilities.hpp index 05bfaf2ac..1a4a69c79 100644 --- a/opm/core/utility/miscUtilities.hpp +++ b/opm/core/utility/miscUtilities.hpp @@ -188,6 +188,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, @@ -320,4 +340,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..9f1e693fe --- /dev/null +++ b/opm/core/utility/miscUtilities_impl.hpp @@ -0,0 +1,41 @@ +#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; + } + } + } + } + } +} diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 2241cacd6..85893f220 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include @@ -182,13 +183,17 @@ namespace { using namespace std; std::array cube; - int num_local_faces = grid.cell_facepos[cell + 1] - grid.cell_facepos[cell]; + typedef Opm::UgGridHelpers::Cell2FacesTraits::Type Cell2Faces; + typedef Cell2Faces::row_type FaceRow; + Cell2Faces c2f=Opm::UgGridHelpers::cell2Faces(grid); + FaceRow faces=c2f[cell]; + int num_local_faces = faces.size(); 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); @@ -417,10 +424,11 @@ 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) { + if (UgGridHelpers::dimensions(grid) != 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. @@ -489,17 +497,17 @@ namespace Opm // 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; + const int* global_cell = UgGridHelpers::globalCell(grid); + const int* cpgdim = UgGridHelpers::cartDims(grid); std::map cartesian_to_compressed; if (global_cell) { - for (int i = 0; i < grid.number_of_cells; ++i) { + for (int i = 0; i < UgGridHelpers::numCells(grid); ++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 < UgGridHelpers::numCells(grid); ++i) { cartesian_to_compressed.insert(std::make_pair(i, i)); } } @@ -553,8 +561,10 @@ namespace Opm 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 << ')'); + 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; @@ -568,7 +578,8 @@ namespace Opm 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]; + int dimensions = UgGridHelpers::dimensions(grid); + const double* cell_perm = &permeability[dimensions*dimensions*cell]; pd.well_index = computeWellIndex(radius, cubical, cell_perm, compdat.compdat[kw].skin_factor_); } @@ -586,7 +597,7 @@ namespace Opm // Set up reference depths that were defaulted. Count perfs. int num_perfs = 0; - assert(grid.dimensions == 3); + assert(UgGridHelpers::dimensions(grid) == 3); for (int w = 0; w < num_wells; ++w) { num_perfs += wellperf_data[w].size(); if (well_data[w].reference_bhp_depth < 0.0) { @@ -594,7 +605,8 @@ namespace Opm 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]; + double depth = UgGridHelpers:: + cellCentroidCoordinate(grid, wellperf_data[w][perf].cell, 2); min_depth = std::min(min_depth, depth); } well_data[w].reference_bhp_depth = min_depth; @@ -1058,18 +1070,18 @@ 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)); } } @@ -1118,12 +1130,14 @@ namespace Opm int j = completion->getJ(); int k = completion->getK(); - const int* cpgdim = grid.cartdims; + const int* cpgdim = UgGridHelpers::cartDims(grid); 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() << ')'); + 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; @@ -1137,7 +1151,8 @@ namespace Opm 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]; + int dimensions=UgGridHelpers::dimensions(grid); + const double* cell_perm = &permeability[dimensions*dimensions*cell]; pd.well_index = computeWellIndex(radius, cubical, cell_perm, completion->getDiameter()); } wellperf_data[well_index].push_back(pd); @@ -1151,7 +1166,7 @@ namespace Opm const int num_wells = well_data.size(); int num_perfs = 0; - assert(grid.dimensions == 3); + assert(=UgGridHelpers::dimensions(grid) == 3); for (int w = 0; w < num_wells; ++w) { num_perfs += wellperf_data[w].size(); if (well_data[w].reference_bhp_depth < 0.0) { @@ -1159,7 +1174,8 @@ namespace Opm 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]; + double depth = UgGridHelpers:: + cellCentroidCoordinate(grid,wellperf_data[w][perf].cell, 2); min_depth = std::min(min_depth, depth); } well_data[w].reference_bhp_depth = min_depth; diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp index 77bbb9038..7a21cd318 100644 --- a/opm/core/wells/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -65,7 +65,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 @@ -73,14 +73,16 @@ 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); WellsManager(const Opm::EclipseStateConstPtr eclipseState, const size_t timeStep, const Opm::EclipseGridParser& deck, const UnstructuredGrid& grid, - const double* permeability); + const double* permeability, + bool checkCellExistence=true); /// Destructor. ~WellsManager(); @@ -135,9 +137,9 @@ namespace Opm private: // 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); @@ -155,6 +157,7 @@ namespace Opm // Data Wells* w_; WellCollection well_collection_; + bool checkCellExistence_; }; } // namespace Opm From b328c13bc0ba2a59402800f253343e8eb022c9b0 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Thu, 20 Feb 2014 10:45:08 +0100 Subject: [PATCH 2/8] Assume begin_cell_centroid to what the name tells us: an iterator over then centroids. Therfore we do not need to call center() in getCoordinate. This is done in the iterator class returned from CpGrid. --- opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp index f180662b0..0715d1741 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp @@ -494,7 +494,7 @@ namespace template double getCoordinate(T t, int i) { - return t->center()[i]; + return (*t)[i]; } From c7e20e5eb3ea8504de71b1607fc1c9645fafed0c Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 25 Feb 2014 15:10:18 +0100 Subject: [PATCH 3/8] Fixed unused parameter warning in initBlackoilStateFromDeck. --- opm/core/simulator/initState.hpp | 1 - opm/core/simulator/initState_impl.hpp | 3 +-- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/opm/core/simulator/initState.hpp b/opm/core/simulator/initState.hpp index 5841592a6..6bdbbef5d 100644 --- a/opm/core/simulator/initState.hpp +++ b/opm/core/simulator/initState.hpp @@ -217,7 +217,6 @@ namespace Opm template void initBlackoilStateFromDeck(int number_of_cells, const int* global_cell, - const int* cartdims, int number_of_faces, FaceCells face_cells, FCI begin_face_centroids, diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index e70c97090..5bcbbd486 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -953,7 +953,7 @@ namespace Opm const double gravity, State& state) { - initBlackoilStateFromDeck(grid.number_of_cells, grid.global_cell, grid.cartdims, + 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); @@ -962,7 +962,6 @@ namespace Opm template void initBlackoilStateFromDeck(int number_of_cells, const int* global_cell, - const int* cartdims, int number_of_faces, FaceCells face_cells, FCI begin_face_centroids, From e30031dc777f1ffe60c3cb46ee04336f8f6da1a0 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 25 Feb 2014 15:10:46 +0100 Subject: [PATCH 4/8] Added missing namespace qualifiers. --- opm/core/simulator/initState_impl.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index 5bcbbd486..32758f0e3 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -441,7 +441,7 @@ 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 = getCoordinate(begin_cell_centroids, dimensions - 1); + 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 { @@ -450,7 +450,7 @@ namespace Opm 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 = getCoordinate(begin_cell_centroids, dimensions - 1); + 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); } From ab0b2a9c6b62fb9faeb52c420848669153b68b17 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 25 Feb 2014 15:12:16 +0100 Subject: [PATCH 5/8] Added support for computePoreVolume for grids apart from UnstructuredGrid. --- opm/core/utility/miscUtilities.cpp | 57 +++-------------- opm/core/utility/miscUtilities.hpp | 43 +++++++++++++ opm/core/utility/miscUtilities_impl.hpp | 82 +++++++++++++++++++++++++ 3 files changed, 132 insertions(+), 50 deletions(-) diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index 73dd738ff..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) @@ -482,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 1a4a69c79..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) @@ -238,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 diff --git a/opm/core/utility/miscUtilities_impl.hpp b/opm/core/utility/miscUtilities_impl.hpp index 9f1e693fe..25bc56111 100644 --- a/opm/core/utility/miscUtilities_impl.hpp +++ b/opm/core/utility/miscUtilities_impl.hpp @@ -1,4 +1,7 @@ #include +#include +#include + namespace Opm { /// @brief Estimates a scalar cell velocity from face fluxes. @@ -38,4 +41,83 @@ namespace Opm } } } + + 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); + } + } + } } From d307df71e7e8935dcdf92774c3ab153c9f6aace1 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 25 Feb 2014 15:45:15 +0100 Subject: [PATCH 6/8] Added suport for grid apart from UG for BlackoilPropertiesFromDeck --- opm/core/props/BlackoilPropertiesFromDeck.hpp | 39 +++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index 78bb74d5f..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(); From f356d8670128699e6bfc11d28a904b31f368902e Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 25 Feb 2014 15:45:50 +0100 Subject: [PATCH 7/8] Remove ambiguous functions increment and getCoordinate from unnamed namespace in favor to those from namespace UgGridHelpers --- .../satfunc/SaturationPropsFromDeck_impl.hpp | 35 +++---------------- 1 file changed, 4 insertions(+), 31 deletions(-) diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp index 0715d1741..542581541 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 @@ -470,36 +471,6 @@ namespace Opm return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]]; } -namespace -{ - - template - const T* increment(T* cc, int i, int dim) - { - return cc+(i*dim); - } - - template - T increment(const T& t, int i, int) - { - return t+i; - } - - template - double getCoordinate(T* cc, int i) - { - return cc[i]; - } - - template - double getCoordinate(T t, int i) - { - return (*t)[i]; - } - - -} - // Initialize saturation scaling parameter template @@ -1002,7 +973,9 @@ namespace if (table[itab][jtab][0] != -1.0) { std::vector& depth = table[0][jtab]; std::vector& val = table[itab][jtab]; - double zc = getCoordinate(increment(begin_cell_centroid, cell, dimensions), + 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); From 2d54a7398b17901a695b2800d1cc1af176c2dafe Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 25 Feb 2014 17:52:51 +0100 Subject: [PATCH 8/8] Refactored WellsManager to partial use without UnstructuredGrid. --- opm/core/wells/WellsManager.cpp | 779 +------------------------ opm/core/wells/WellsManager.hpp | 32 +- opm/core/wells/WellsManager_impl.hpp | 827 +++++++++++++++++++++++++++ 3 files changed, 879 insertions(+), 759 deletions(-) create mode 100644 opm/core/wells/WellsManager_impl.hpp diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 47ec79323..3788f3611 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -23,17 +23,14 @@ #include #include #include -#include #include #include #include -#include #include #include #include -#include #include #include #include @@ -45,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() { @@ -123,8 +116,6 @@ namespace namespace InjectionControl { - enum Mode { RATE, RESV, BHP, - THP, GRUP }; namespace Details { std::map @@ -178,32 +169,6 @@ namespace } // namespace InjectionControl - - std::array getCubeDim(const UnstructuredGrid& grid, int cell) - { - using namespace std; - std::array cube; - typedef Opm::UgGridHelpers::Cell2FacesTraits::Type Cell2Faces; - typedef Cell2Faces::row_type FaceRow; - Cell2Faces c2f=Opm::UgGridHelpers::cell2Faces(grid); - FaceRow faces=c2f[cell]; - int num_local_faces = faces.size(); - vector x(num_local_faces); - vector y(num_local_faces); - vector z(num_local_faces); - for (int lf=0; lf 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 = UgGridHelpers::globalCell(grid); - const int* cpgdim = UgGridHelpers::cartDims(grid); - std::map cartesian_to_compressed; - - if (global_cell) { - for (int i = 0; i < UgGridHelpers::numCells(grid); ++i) { - cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); - } - } - else { - for (int i = 0; i < UgGridHelpers::numCells(grid); ++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 = getCubeDim(grid, cell); - int dimensions = UgGridHelpers::dimensions(grid); - const double* cell_perm = &permeability[dimensions*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(UgGridHelpers::dimensions(grid) == 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:: - cellCentroidCoordinate(grid, 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") { - int cur_ctrl = well_controls_get_current(w_->ctrls[index]); - if (cur_ctrl >= 0) { - well_controls_invert_current(w_->ctrls[index]); - } - assert(well_controls_get_current(w_->ctrls[index]) < 0); - } else if (line.openshutflag_ == "OPEN") { - int cur_ctrl = well_controls_get_current(w_->ctrls[index]); - if (cur_ctrl < 0) { - well_controls_invert_current(w_->ctrls[index]); - } - assert(well_controls_get_current(w_->ctrls[index]) >= 0); - } 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() { @@ -1021,127 +408,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) } - 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 = UgGridHelpers::cartDims(grid); - 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 = getCubeDim(grid, cell); - int dimensions=UgGridHelpers::dimensions(grid); - const double* cell_perm = &permeability[dimensions*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(=UgGridHelpers::dimensions(grid) == 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:: - cellCentroidCoordinate(grid,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) { @@ -1159,7 +426,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) 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); @@ -1179,7 +446,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) } 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); @@ -1199,7 +466,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) } 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, @@ -1218,9 +485,9 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) { - 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]); } @@ -1266,7 +533,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) 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, @@ -1280,7 +547,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) 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, @@ -1294,7 +561,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) 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, @@ -1311,7 +578,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) 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; @@ -1323,7 +590,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) } 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), @@ -1333,7 +600,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) } 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, @@ -1349,9 +616,9 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence) 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 450812c0a..32a8fb441 100644 --- a/opm/core/wells/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -78,6 +78,18 @@ namespace Opm 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, @@ -136,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, bool checkCellExistence=true); WellsManager& operator=(const WellsManager& other); 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, @@ -164,5 +190,5 @@ namespace Opm } // 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..11de0cb08 --- /dev/null +++ b/opm/core/wells/WellsManager_impl.hpp @@ -0,0 +1,827 @@ +#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") { + int cur_ctrl = well_controls_get_current(w_->ctrls[index]); + if (cur_ctrl >= 0) { + well_controls_invert_current(w_->ctrls[index]); + } + assert(well_controls_get_current(w_->ctrls[index]) < 0); + } else if (line.openshutflag_ == "OPEN") { + int cur_ctrl = well_controls_get_current(w_->ctrls[index]); + if (cur_ctrl < 0) { + well_controls_invert_current(w_->ctrls[index]); + } + assert(well_controls_get_current(w_->ctrls[index]) >= 0); + } 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