From d5f470cb686d08b4080022ae8b12be61708e154d Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Mon, 17 Feb 2014 13:23:01 +0100 Subject: [PATCH] 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