diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index 77aaf4ca7..c32ad02ee 100644 --- a/opm/core/linalg/LinearSolverIstl.cpp +++ b/opm/core/linalg/LinearSolverIstl.cpp @@ -30,6 +30,7 @@ // TODO: clean up includes. #include +#include #include #include #include @@ -39,6 +40,10 @@ #include #include +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) +#include +#endif + #include #include @@ -61,7 +66,7 @@ namespace Opm solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity, double prolongateFactor, int smoothsteps); -#ifdef HAS_DUNE_FAST_AMG +#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) LinearSolverInterface::LinearSolverReport solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity, double prolongateFactor, int smoothsteps); @@ -173,7 +178,7 @@ namespace Opm linsolver_prolongate_factor_, linsolver_smooth_steps_); break; case KAMG: -#ifdef HAS_DUNE_FAST_AMG +#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) res = solveKAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_, linsolver_prolongate_factor_, linsolver_smooth_steps_); #else @@ -181,7 +186,7 @@ namespace Opm #endif break; case FastAMG: -#ifdef HAS_DUNE_FAST_AMG +#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) res = solveFastAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_, linsolver_prolongate_factor_); #else @@ -311,7 +316,7 @@ namespace Opm } -#ifdef HAS_DUNE_FAST_AMG +#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) LinearSolverInterface::LinearSolverReport solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity, double linsolver_prolongate_factor, int linsolver_smooth_steps) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index 4c2ca10ed..b773408f5 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -869,7 +869,7 @@ assemble_completion_to_well(int i, int w, int c, int nc, int np, W = wells->W; ctrl = W->ctrls[ w ]; - if (well_controls_get_current(ctrl) < 0) { + if (well_controls_well_is_shut( ctrl )) { /* Interpreting a negative current control index to mean a shut well */ welleq_coeff_shut(np, h, &res, &w2c, &w2w); } @@ -933,7 +933,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells , for (w = i = 0; w < W->number_of_wells; w++) { pw = wpress[ w ]; - is_open = (well_controls_get_current(W->ctrls[w]) >= 0); + is_open = well_controls_well_is_open(W->ctrls[w]); for (; i < W->well_connpos[w + 1]; i++, pmobp += np) { diff --git a/opm/core/pressure/tpfa/ifs_tpfa.c b/opm/core/pressure/tpfa/ifs_tpfa.c index 5130ba0b3..a9e2ce73c 100644 --- a/opm/core/pressure/tpfa/ifs_tpfa.c +++ b/opm/core/pressure/tpfa/ifs_tpfa.c @@ -363,7 +363,7 @@ assemble_well_contrib(int nc , for (w = 0; w < W->number_of_wells; w++) { ctrls = W->ctrls[ w ]; - if (well_controls_get_current(ctrls) < 0) { + if (well_controls_well_is_shut(ctrls) ) { /* Treat this well as a shut well, isolated from the domain. */ diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp index caa27f1b9..a7de41c8b 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp @@ -23,7 +23,6 @@ namespace Opm { - BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, const UnstructuredGrid& grid, bool init_rock) @@ -43,6 +42,25 @@ namespace Opm } } + BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + bool init_rock) + { + if (init_rock){ + rock_.init(newParserDeck, grid); + } + pvt_.init(newParserDeck, /*numSamples=*/0); + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, /*numSamples=*/0); + + if (pvt_.numPhases() != satprops_->numPhases()) { + OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); + } + } + BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, const UnstructuredGrid& grid, const parameter::ParameterGroup& param, @@ -58,8 +76,8 @@ namespace Opm // Unfortunate lack of pointer smartness here... const int sat_samples = param.getDefault("sat_tab_size", 0); std::string threephase_model = param.getDefault("threephase_model", "simple"); - if (deck.hasField("ENDSCALE") && threephase_model != "simple") { - OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'simple' model only."); + if (deck.hasField("ENDSCALE") && threephase_model != "gwseg") { + OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'gwseg' model only."); } if (sat_samples > 1) { if (threephase_model == "stone2") { @@ -107,6 +125,70 @@ namespace Opm } } + BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const parameter::ParameterGroup& param, + bool init_rock) + { + if(init_rock){ + rock_.init(newParserDeck, grid); + } + + const int pvt_samples = param.getDefault("pvt_tab_size", 0); + pvt_.init(newParserDeck, pvt_samples); + + // Unfortunate lack of pointer smartness here... + const int sat_samples = param.getDefault("sat_tab_size", 0); + std::string threephase_model = param.getDefault("threephase_model", "gwseg"); + if (newParserDeck->hasKeyword("ENDSCALE") && threephase_model != "gwseg") { + OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'gwseg' model only."); + } + if (sat_samples > 1) { + if (threephase_model == "stone2") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else if (threephase_model == "simple") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else if (threephase_model == "gwseg") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else { + OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model); + } + } else { + if (threephase_model == "stone2") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else if (threephase_model == "simple") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else if (threephase_model == "gwseg") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else { + OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model); + } + } + + if (pvt_.numPhases() != satprops_->numPhases()) { + OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); + } + } + BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck() { } diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index 49072b29e..5c74a76a6 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -27,6 +27,9 @@ #include #include #include + +#include + #include struct UnstructuredGrid; @@ -44,7 +47,15 @@ namespace Opm /// \param[in] grid Grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. - BlackoilPropertiesFromDeck(const EclipseGridParser& deck, + BlackoilPropertiesFromDeck(const EclipseGridParser &deck, + const UnstructuredGrid& grid, bool init_rock=true ); + + /// Initialize from deck and grid. + /// \param[in] deck Deck input parser + /// \param[in] grid Grid to which property object applies, needed for the + /// mapping from cell indices (typically from a processed grid) + /// to logical cartesian indices consistent with the deck. + BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, const UnstructuredGrid& grid, bool init_rock=true ); /// Initialize from deck, grid and parameters. @@ -63,6 +74,22 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock=true); + /// Initialize from deck, grid and parameters. + /// \param[in] deck Deck input parser + /// \param[in] grid Grid to which property object applies, needed for the + /// mapping from cell indices (typically from a processed grid) + /// to logical cartesian indices consistent with the deck. + /// \param[in] param Parameters. Accepted parameters include: + /// pvt_tab_size (200) number of uniform sample points for dead-oil pvt tables. + /// sat_tab_size (200) number of uniform sample points for saturation tables. + /// threephase_model("simple") three-phase relperm model (accepts "simple" and "stone2"). + /// For both size parameters, a 0 or negative value indicates that no spline fitting is to + /// be done, and the input fluid data used directly for linear interpolation. + BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const parameter::ParameterGroup& param, + bool init_rock=true); + /// Destructor. virtual ~BlackoilPropertiesFromDeck(); diff --git a/opm/core/props/phaseUsageFromDeck.hpp b/opm/core/props/phaseUsageFromDeck.hpp index 01f28ee14..6ee66075e 100644 --- a/opm/core/props/phaseUsageFromDeck.hpp +++ b/opm/core/props/phaseUsageFromDeck.hpp @@ -24,10 +24,49 @@ #include #include +#include +#include + namespace Opm { + /// Looks at presence of WATER, OIL and GAS keywords in state object + /// to determine active phases. + inline PhaseUsage phaseUsageFromDeck(Opm::EclipseStateConstPtr eclipseState) + { + PhaseUsage pu; + std::fill(pu.phase_used, pu.phase_used + BlackoilPhases::MaxNumPhases, 0); + + // Discover phase usage. + if (eclipseState->hasPhase(Phase::PhaseEnum::WATER)) { + pu.phase_used[BlackoilPhases::Aqua] = 1; + } + if (eclipseState->hasPhase(Phase::PhaseEnum::OIL)) { + pu.phase_used[BlackoilPhases::Liquid] = 1; + } + if (eclipseState->hasPhase(Phase::PhaseEnum::GAS)) { + pu.phase_used[BlackoilPhases::Vapour] = 1; + } + pu.num_phases = 0; + for (int i = 0; i < BlackoilPhases::MaxNumPhases; ++i) { + pu.phase_pos[i] = pu.num_phases; + pu.num_phases += pu.phase_used[i]; + } + + // Only 2 or 3 phase systems handled. + if (pu.num_phases < 2 || pu.num_phases > 3) { + OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases."); + } + + // We need oil systems, since we do not support the keywords needed for + // water-gas systems. + if (!pu.phase_used[BlackoilPhases::Liquid]) { + OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems."); + } + + return pu; + } /// Looks at presence of WATER, OIL and GAS keywords in deck /// to determine active phases. @@ -66,6 +105,43 @@ namespace Opm return pu; } + /// Looks at presence of WATER, OIL and GAS keywords in deck + /// to determine active phases. + inline PhaseUsage phaseUsageFromDeck(Opm::DeckConstPtr newParserDeck) + { + PhaseUsage pu; + std::fill(pu.phase_used, pu.phase_used + BlackoilPhases::MaxNumPhases, 0); + + // Discover phase usage. + if (newParserDeck->hasKeyword("WATER")) { + pu.phase_used[BlackoilPhases::Aqua] = 1; + } + if (newParserDeck->hasKeyword("OIL")) { + pu.phase_used[BlackoilPhases::Liquid] = 1; + } + if (newParserDeck->hasKeyword("GAS")) { + pu.phase_used[BlackoilPhases::Vapour] = 1; + } + pu.num_phases = 0; + for (int i = 0; i < BlackoilPhases::MaxNumPhases; ++i) { + pu.phase_pos[i] = pu.num_phases; + pu.num_phases += pu.phase_used[i]; + } + + // Only 2 or 3 phase systems handled. + if (pu.num_phases < 2 || pu.num_phases > 3) { + OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases."); + } + + // We need oil systems, since we do not support the keywords needed for + // water-gas systems. + if (!pu.phase_used[BlackoilPhases::Liquid]) { + OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems."); + } + + return pu; + } + } #endif // OPM_PHASEUSAGEFROMDECK_HEADER_INCLUDED diff --git a/opm/core/props/rock/RockCompressibility.cpp b/opm/core/props/rock/RockCompressibility.cpp index 7b0f69cb8..d1729533c 100644 --- a/opm/core/props/rock/RockCompressibility.cpp +++ b/opm/core/props/rock/RockCompressibility.cpp @@ -25,6 +25,9 @@ #include #include +#include +#include + #include namespace Opm @@ -65,6 +68,44 @@ namespace Opm } } + RockCompressibility::RockCompressibility(Opm::DeckConstPtr newParserDeck) + : pref_(0.0), + rock_comp_(0.0) + { + if (newParserDeck->hasKeyword("ROCKTAB")) { + Opm::DeckKeywordConstPtr rtKeyword = newParserDeck->getKeyword("ROCKTAB"); + if (rtKeyword->size() != 1) + OPM_THROW(std::runtime_error, "Can only handle a single region in ROCKTAB."); + + // the number of colums of the "ROCKTAB" keyword + // depends on the presence of the "RKTRMDIR" + // keyword. Messy stuff... + bool isDirectional = newParserDeck->hasKeyword("RKTRMDIR"); + if (isDirectional) + { + // well, okay. we don't support non-isotropic + // transmissibility multipliers yet + OPM_THROW(std::runtime_error, "Support for non-isotropic " + "transmissibility multipliers is not implemented yet."); + }; + + Opm::RocktabTable rocktabTable(rtKeyword, isDirectional); + + p_ = rocktabTable.getPressureColumn(); + poromult_ = rocktabTable.getPoreVolumeMultiplierColumn(); + transmult_ = rocktabTable.getTransmissibilityMultiplierColumn(); + } else if (newParserDeck->hasKeyword("ROCK")) { + Opm::RockTable rockTable(newParserDeck->getKeyword("ROCK")); + if (rockTable.numRows() != 1) + OPM_THROW(std::runtime_error, "Can only handle a single region in ROCK."); + + pref_ = rockTable.getPressureColumn()[0]; + rock_comp_ = rockTable.getCompressibilityColumn()[0]; + } else { + std::cout << "**** warning: no rock compressibility data found in deck (ROCK or ROCKTAB)." << std::endl; + } + } + bool RockCompressibility::isActive() const { return !p_.empty() || (rock_comp_ != 0.0); diff --git a/opm/core/props/rock/RockCompressibility.hpp b/opm/core/props/rock/RockCompressibility.hpp index 03567caed..bd375e1ca 100644 --- a/opm/core/props/rock/RockCompressibility.hpp +++ b/opm/core/props/rock/RockCompressibility.hpp @@ -20,6 +20,8 @@ #ifndef OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED #define OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED +#include + #include namespace Opm @@ -35,6 +37,10 @@ namespace Opm /// Looks for the keywords ROCK and ROCKTAB. RockCompressibility(const EclipseGridParser& deck); + /// Construct from input deck. + /// Looks for the keywords ROCK and ROCKTAB. + RockCompressibility(Opm::DeckConstPtr newParserDeck); + /// Construct from parameters. /// Accepts the following parameters (with defaults). /// rock_compressibility_pref (100.0) [given in bar] diff --git a/opm/core/props/rock/RockFromDeck.cpp b/opm/core/props/rock/RockFromDeck.cpp index 730363434..1058fa997 100644 --- a/opm/core/props/rock/RockFromDeck.cpp +++ b/opm/core/props/rock/RockFromDeck.cpp @@ -21,6 +21,9 @@ #include "config.h" #include #include + +#include + #include namespace Opm @@ -37,6 +40,10 @@ namespace Opm PermeabilityKind fillTensor(const EclipseGridParser& parser, std::vector*>& tensor, std::array& kmap); + PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck, + std::vector*>& tensor, + std::array& kmap); + } // anonymous namespace @@ -64,6 +71,15 @@ namespace Opm assignPermeability(deck, grid, perm_threshold); } + void RockFromDeck::init(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid) + { + assignPorosity(newParserDeck, grid); + permfield_valid_.assign(grid.number_of_cells, false); + const double perm_threshold = 0.0; // Maybe turn into parameter? + assignPermeability(newParserDeck, grid, perm_threshold); + } + void RockFromDeck::assignPorosity(const EclipseGridParser& parser, const UnstructuredGrid& grid) @@ -79,6 +95,21 @@ namespace Opm } } + void RockFromDeck::assignPorosity(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid) + { + porosity_.assign(grid.number_of_cells, 1.0); + const int* gc = grid.global_cell; + if (newParserDeck->hasKeyword("PORO")) { + const std::vector& poro = newParserDeck->getKeyword("PORO")->getSIDoubleData(); + for (int c = 0; c < int(porosity_.size()); ++c) { + const int deck_pos = (gc == NULL) ? c : gc[c]; + assert(0 <= c && c < (int) porosity_.size()); + assert(0 <= deck_pos && deck_pos < (int) poro.size()); + porosity_[c] = poro[deck_pos]; + } + } + } void RockFromDeck::assignPermeability(const EclipseGridParser& parser, @@ -135,6 +166,59 @@ namespace Opm } } + void RockFromDeck::assignPermeability(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + 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; + + assert(num_global_cells > 0); + + permeability_.assign(dim * dim * nc, 0.0); + + std::vector*> tensor; + tensor.reserve(10); + + const std::vector zero(num_global_cells, 0.0); + tensor.push_back(&zero); + + std::array kmap; + PermeabilityKind pkind = fillTensor(newParserDeck, tensor, kmap); + if (pkind == Invalid) { + OPM_THROW(std::runtime_error, "Invalid permeability field."); + } + + // Assign permeability values only if such values are + // given in the input deck represented by 'newParserDeck'. In + // other words: Don't set any (arbitrary) default values. + // It is infinitely better to experience a reproducible + // crash than subtle errors resulting from a (poorly + // 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) { + // SharedPermTensor K(dim, dim, &permeability_[off]); + int kix = 0; + const int glob = (gc == NULL) ? c : gc[c]; + + for (int i = 0; i < dim; ++i) { + for (int j = 0; j < dim; ++j, ++kix) { + // K(i,j) = (*tensor[kmap[kix]])[glob]; + permeability_[off + kix] = (*tensor[kmap[kix]])[glob]; + } + // K(i,i) = std::max(K(i,i), perm_threshold); + permeability_[off + 3*i + i] = std::max(permeability_[off + 3*i + i], perm_threshold); + } + + permfield_valid_[c] = std::vector::value_type(1); + } + } + } namespace { @@ -204,6 +288,72 @@ namespace Opm return retval; } + /// @brief + /// Classify and verify a given permeability specification + /// from a structural point of view. In particular, we + /// verify that there are no off-diagonal permeability + /// components such as @f$k_{xy}@f$ unless the + /// corresponding diagonal components are known as well. + /// + /// @param newParserDeck [in] + /// An Eclipse data parser capable of answering which + /// permeability components are present in a given input + /// deck. + /// + /// @return + /// An enum value with the following possible values: + /// ScalarPerm only one component was given. + /// DiagonalPerm more than one component given. + /// TensorPerm at least one cross-component given. + /// None no components given. + /// Invalid invalid set of components given. + PermeabilityKind classifyPermeability(Opm::DeckConstPtr newParserDeck) + { + const bool xx = newParserDeck->hasKeyword("PERMX" ); + const bool xy = newParserDeck->hasKeyword("PERMXY"); + const bool xz = newParserDeck->hasKeyword("PERMXZ"); + + const bool yx = newParserDeck->hasKeyword("PERMYX"); + const bool yy = newParserDeck->hasKeyword("PERMY" ); + const bool yz = newParserDeck->hasKeyword("PERMYZ"); + + const bool zx = newParserDeck->hasKeyword("PERMZX"); + const bool zy = newParserDeck->hasKeyword("PERMZY"); + const bool zz = newParserDeck->hasKeyword("PERMZ" ); + + int num_cross_comp = xy + xz + yx + yz + zx + zy; + int num_comp = xx + yy + zz + num_cross_comp; + PermeabilityKind retval = None; + if (num_cross_comp > 0) { + retval = TensorPerm; + } else { + if (num_comp == 1) { + retval = ScalarPerm; + } else if (num_comp >= 2) { + retval = DiagonalPerm; + } + } + + bool ok = true; + if (num_comp > 0) { + // At least one tensor component specified on input. + // Verify that any remaining components are OK from a + // structural point of view. In particular, there + // must not be any cross-components (e.g., k_{xy}) + // unless the corresponding diagonal component (e.g., + // k_{xx}) is present as well... + // + ok = xx || !(xy || xz || yx || zx) ; + ok = ok && (yy || !(yx || yz || xy || zy)); + ok = ok && (zz || !(zx || zy || xz || yz)); + } + if (!ok) { + retval = Invalid; + } + + return retval; + } + /// @brief /// Copy isotropic (scalar) permeability to other diagonal @@ -333,6 +483,106 @@ namespace Opm return kind; } + /// @brief + /// Extract pointers to appropriate tensor components from + /// input deck. The permeability tensor is, generally, + /// @code + /// [ kxx kxy kxz ] + /// K = [ kyx kyy kyz ] + /// [ kzx kzy kzz ] + /// @endcode + /// We store these values in a linear array using natural + /// ordering with the column index cycling the most rapidly. + /// In particular we use the representation + /// @code + /// [ 0 1 2 3 4 5 6 7 8 ] + /// K = [ kxx, kxy, kxz, kyx, kyy, kyz, kzx, kzy, kzz ] + /// @endcode + /// Moreover, we explicitly enforce symmetric tensors by + /// assigning + /// @code + /// 3 1 6 2 7 5 + /// kyx = kxy, kzx = kxz, kzy = kyz + /// @endcode + /// However, we make no attempt at enforcing positive + /// definite tensors. + /// + /// @param [in] parser + /// An Eclipse data parser capable of answering which + /// permeability components are present in a given input + /// deck as well as retrieving the numerical value of each + /// permeability component in each grid cell. + /// + /// @param [out] tensor + /// @param [out] kmap + PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck, + std::vector*>& tensor, + std::array& kmap) + { + PermeabilityKind kind = classifyPermeability(newParserDeck); + if (kind == Invalid) { + OPM_THROW(std::runtime_error, "Invalid set of permeability fields given."); + } + assert(tensor.size() == 1); + for (int i = 0; i < 9; ++i) { kmap[i] = 0; } + + enum { xx, xy, xz, // 0, 1, 2 + yx, yy, yz, // 3, 4, 5 + zx, zy, zz }; // 6, 7, 8 + + // ----------------------------------------------------------- + // 1st row: [kxx, kxy, kxz] + if (newParserDeck->hasKeyword("PERMX" )) { + kmap[xx] = tensor.size(); + tensor.push_back(&newParserDeck->getKeyword("PERMX")->getSIDoubleData()); + + setScalarPermIfNeeded(kmap, xx, yy, zz); + } + if (newParserDeck->hasKeyword("PERMXY")) { + kmap[xy] = kmap[yx] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMXY")->getSIDoubleData()); + } + if (newParserDeck->hasKeyword("PERMXZ")) { + kmap[xz] = kmap[zx] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMXZ")->getSIDoubleData()); + } + + // ----------------------------------------------------------- + // 2nd row: [kyx, kyy, kyz] + if (newParserDeck->hasKeyword("PERMYX")) { + kmap[yx] = kmap[xy] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMYX")->getSIDoubleData()); + } + if (newParserDeck->hasKeyword("PERMY" )) { + kmap[yy] = tensor.size(); + tensor.push_back(&newParserDeck->getKeyword("PERMY")->getSIDoubleData()); + + setScalarPermIfNeeded(kmap, yy, zz, xx); + } + if (newParserDeck->hasKeyword("PERMYZ")) { + kmap[yz] = kmap[zy] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMYZ")->getSIDoubleData()); + } + + // ----------------------------------------------------------- + // 3rd row: [kzx, kzy, kzz] + if (newParserDeck->hasKeyword("PERMZX")) { + kmap[zx] = kmap[xz] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMZX")->getSIDoubleData()); + } + if (newParserDeck->hasKeyword("PERMZY")) { + kmap[zy] = kmap[yz] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMZY")->getSIDoubleData()); + } + if (newParserDeck->hasKeyword("PERMZ" )) { + kmap[zz] = tensor.size(); + tensor.push_back(&newParserDeck->getKeyword("PERMZ")->getSIDoubleData()); + + setScalarPermIfNeeded(kmap, zz, xx, yy); + } + return kind; + } + } // anonymous namespace } // namespace Opm diff --git a/opm/core/props/rock/RockFromDeck.hpp b/opm/core/props/rock/RockFromDeck.hpp index e277ee2a5..b5e0fb4e0 100644 --- a/opm/core/props/rock/RockFromDeck.hpp +++ b/opm/core/props/rock/RockFromDeck.hpp @@ -22,6 +22,9 @@ #include + +#include + #include struct UnstructuredGrid; @@ -43,6 +46,14 @@ namespace Opm void init(const EclipseGridParser& deck, const UnstructuredGrid& grid); + /// Initialize from deck and grid. + /// \param newParserDeck Deck produced by the opm-parser code + /// \param grid Grid to which property object applies, needed for the + /// mapping from cell indices (typically from a processed grid) + /// to logical cartesian indices consistent with the deck. + void init(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid); + /// \return D, the number of spatial dimensions. Always 3 for deck input. int numDimensions() const { @@ -72,9 +83,14 @@ namespace Opm private: void assignPorosity(const EclipseGridParser& parser, const UnstructuredGrid& grid); + void assignPorosity(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid); void assignPermeability(const EclipseGridParser& parser, const UnstructuredGrid& grid, const double perm_threshold); + void assignPermeability(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + double perm_threshold); std::vector porosity_; std::vector permeability_; diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp index c321a42af..757b1481a 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp @@ -27,6 +27,9 @@ #include #include #include + +#include + #include struct UnstructuredGrid; @@ -60,6 +63,17 @@ namespace Opm const UnstructuredGrid& grid, const int samples); + /// Initialize from deck and grid. + /// \param[in] deck Deck input parser + /// \param[in] grid Grid to which property object applies, needed for the + /// mapping from cell indices (typically from a processed grid) + /// to logical cartesian indices consistent with the deck. + /// \param[in] samples Number of uniform sample points for saturation tables. + /// NOTE: samples will only be used with the SatFuncSetUniform template argument. + void init(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const int samples); + /// \return P, the number of phases. int numPhases() const; @@ -132,6 +146,14 @@ namespace Opm const UnstructuredGrid& grid, const std::string& keyword, std::vector& scaleparam); + void initEPS(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid); + void initEPSHyst(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid); + void initEPSKey(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const std::string& keyword, + std::vector& scaleparam); void initEPSParam(const int cell, EPSTransforms::Transform& data, const bool oil, @@ -149,6 +171,11 @@ namespace Opm const std::vector& s0, const std::vector& krsr, const std::vector& krmax); + + bool columnIsMasked_(Opm::DeckConstPtr newParserDeck, + const std::string& keywordName, + int /* columnIdx */) + { return newParserDeck->getKeyword(keywordName)->getRecord(0)->getItem(0)->getSIDouble(0) != -1.0; } }; diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp index c30069934..a6e699faf 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp @@ -26,6 +26,11 @@ #include #include +#include +#include +#include +#include + #include namespace Opm @@ -140,6 +145,123 @@ namespace Opm } + /// Initialize from deck. + template + void SaturationPropsFromDeck::init(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const int samples) + { + phase_usage_ = phaseUsageFromDeck(newParserDeck); + + // Extract input data. + // Oil phase should be active. + if (!phase_usage_.phase_used[Liquid]) { + OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- oil phase must be active."); + } + + // Obtain SATNUM, if it exists, and create cell_to_func_. + // Otherwise, let the cell_to_func_ mapping be just empty. + int satfuncs_expected = 1; + if (newParserDeck->hasKeyword("SATNUM")) { + const std::vector& satnum = newParserDeck->getKeyword("SATNUM")->getIntData(); + satfuncs_expected = *std::max_element(satnum.begin(), satnum.end()); + const int num_cells = grid.number_of_cells; + cell_to_func_.resize(num_cells); + const int* gc = grid.global_cell; + for (int cell = 0; cell < num_cells; ++cell) { + const int newParserDeck_pos = (gc == NULL) ? cell : gc[cell]; + cell_to_func_[cell] = satnum[newParserDeck_pos] - 1; + } + } + + // Find number of tables, check for consistency. + enum { Uninitialized = -1 }; + int num_tables = Uninitialized; + if (phase_usage_.phase_used[Aqua]) { + num_tables = newParserDeck->getKeyword("SWOF")->size(); + if (num_tables < satfuncs_expected) { + OPM_THROW(std::runtime_error, "Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected); + } + } + if (phase_usage_.phase_used[Vapour]) { + int num_sgof_tables = newParserDeck->getKeyword("SGOF")->size(); + if (num_sgof_tables < satfuncs_expected) { + OPM_THROW(std::runtime_error, "Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected); + } + if (num_tables == Uninitialized) { + num_tables = num_sgof_tables; + } else if (num_tables != num_sgof_tables) { + OPM_THROW(std::runtime_error, "Inconsistent number of tables in SWOF and SGOF."); + } + } + + // Initialize tables. + satfuncset_.resize(num_tables); + for (int table = 0; table < num_tables; ++table) { + satfuncset_[table].init(newParserDeck, table, phase_usage_, samples); + } + + // Saturation table scaling + do_hyst_ = false; + do_eps_ = false; + do_3pt_ = false; + if (newParserDeck->hasKeyword("ENDSCALE")) { + Opm::EndscaleWrapper endscale(newParserDeck->getKeyword("ENDSCALE")); + if (endscale.directionSwitch() != std::string("NODIR")) { + OPM_THROW(std::runtime_error, + "SaturationPropsFromDeck::init() -- ENDSCALE: " + "Currently only 'NODIR' accepted."); + } + if (!endscale.isReversible()) { + OPM_THROW(std::runtime_error, + "SaturationPropsFromDeck::init() -- ENDSCALE: " + "Currently only 'REVERS' accepted."); + } + if (newParserDeck->hasKeyword("SCALECRS")) { + Opm::ScalecrsWrapper scalecrs(newParserDeck->getKeyword("SCALECRS")); + if (scalecrs.isEnabled()) { + do_3pt_ = true; + } + } + do_eps_ = true; + + initEPS(newParserDeck, grid); + + // For now, a primitive detection of hysteresis. TODO: SATOPTS HYSTER/ and EHYSTR + do_hyst_ = + newParserDeck->hasKeyword("ISWL") + || newParserDeck->hasKeyword("ISWU") + || newParserDeck->hasKeyword("ISWCR") + || newParserDeck->hasKeyword("ISGL") + || newParserDeck->hasKeyword("ISGU") + || newParserDeck->hasKeyword("ISGCR") + || newParserDeck->hasKeyword("ISOWCR") + || newParserDeck->hasKeyword("ISOGCR"); + if (do_hyst_) { + if (newParserDeck->hasKeyword("KRW") + || newParserDeck->hasKeyword("KRG") + || newParserDeck->hasKeyword("KRO") + || newParserDeck->hasKeyword("KRWR") + || newParserDeck->hasKeyword("KRGR") + || newParserDeck->hasKeyword("KRORW") + || newParserDeck->hasKeyword("KRORG") + || newParserDeck->hasKeyword("IKRW") + || newParserDeck->hasKeyword("IKRG") + || newParserDeck->hasKeyword("IKRO") + || newParserDeck->hasKeyword("IKRWR") + || newParserDeck->hasKeyword("IKRGR") + || newParserDeck->hasKeyword("IKRORW") + || newParserDeck->hasKeyword("IKRORG") ) { + OPM_THROW(std::runtime_error, + "SaturationPropsFromDeck::init() -- ENDSCALE: " + "Currently hysteresis and relperm value scaling " + "cannot be combined."); + } + initEPSHyst(newParserDeck, grid); + } + } + } + /// \return P, the number of phases. @@ -383,6 +505,69 @@ namespace Opm } } + // Initialize saturation scaling parameters + template + void SaturationPropsFromDeck::initEPS(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid) + { + std::vector swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr; + std::vector krw, krg, kro, krwr, krgr, krorw, krorg; + // Initialize saturation scaling parameter + initEPSKey(newParserDeck, grid, std::string("SWL"), swl); + initEPSKey(newParserDeck, grid, std::string("SWU"), swu); + initEPSKey(newParserDeck, grid, std::string("SWCR"), swcr); + initEPSKey(newParserDeck, grid, std::string("SGL"), sgl); + initEPSKey(newParserDeck, grid, std::string("SGU"), sgu); + initEPSKey(newParserDeck, grid, std::string("SGCR"), sgcr); + initEPSKey(newParserDeck, grid, std::string("SOWCR"), sowcr); + initEPSKey(newParserDeck, grid, std::string("SOGCR"), sogcr); + initEPSKey(newParserDeck, grid, std::string("KRW"), krw); + initEPSKey(newParserDeck, grid, std::string("KRG"), krg); + initEPSKey(newParserDeck, grid, std::string("KRO"), kro); + initEPSKey(newParserDeck, grid, std::string("KRWR"), krwr); + initEPSKey(newParserDeck, grid, std::string("KRGR"), krgr); + initEPSKey(newParserDeck, grid, std::string("KRORW"), krorw); + initEPSKey(newParserDeck, grid, std::string("KRORG"), krorg); + + eps_transf_.resize(grid.number_of_cells); + + const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour]; + const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + + for (int cell = 0; cell < grid.number_of_cells; ++cell) { + if (oilWater) { + // ### krw + initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], + funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw); + // ### krow + initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro); + } else if (oilGas) { + // ### krg + initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], + funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg); + // ### krog + initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro); + } else if (threephase) { + // ### krw + initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_, + funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw); + // ### krow + initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, + funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro); + // ### krg + initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_, + funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg); + // ### krog + initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, + funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro); + } + } + } // Initialize hysteresis saturation scaling parameters template @@ -450,6 +635,71 @@ namespace Opm } } + // Initialize hysteresis saturation scaling parameters + template + void SaturationPropsFromDeck::initEPSHyst(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid) + { + std::vector iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr; + std::vector ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg; + // Initialize hysteresis saturation scaling parameters + initEPSKey(newParserDeck, grid, std::string("ISWL"), iswl); + initEPSKey(newParserDeck, grid, std::string("ISWU"), iswu); + initEPSKey(newParserDeck, grid, std::string("ISWCR"), iswcr); + initEPSKey(newParserDeck, grid, std::string("ISGL"), isgl); + initEPSKey(newParserDeck, grid, std::string("ISGU"), isgu); + initEPSKey(newParserDeck, grid, std::string("ISGCR"), isgcr); + initEPSKey(newParserDeck, grid, std::string("ISOWCR"), isowcr); + initEPSKey(newParserDeck, grid, std::string("ISOGCR"), isogcr); + initEPSKey(newParserDeck, grid, std::string("IKRW"), ikrw); + initEPSKey(newParserDeck, grid, std::string("IKRG"), ikrg); + initEPSKey(newParserDeck, grid, std::string("IKRO"), ikro); + initEPSKey(newParserDeck, grid, std::string("IKRWR"), ikrwr); + initEPSKey(newParserDeck, grid, std::string("IKRGR"), ikrgr); + initEPSKey(newParserDeck, grid, std::string("IKRORW"), ikrorw); + initEPSKey(newParserDeck, grid, std::string("IKRORG"), ikrorg); + + eps_transf_hyst_.resize(grid.number_of_cells); + sat_hyst_.resize(grid.number_of_cells); + + const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour]; + const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + + for (int cell = 0; cell < grid.number_of_cells; ++cell) { + if (oilWater) { + // ### krw + initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], + funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw); + // ### krow + initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro); + } else if (oilGas) { + // ### krg + initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], + funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg); + // ### krog + initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro); + } else if (threephase) { + // ### krw + initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_, + funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw); + // ### krow + initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, + funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro); + // ### krg + initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_, + funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg); + // ### krog + initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, + funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro); + } + } + } + // Initialize saturation scaling parameter template void SaturationPropsFromDeck::initEPSKey(const EclipseGridParser& deck, @@ -624,6 +874,213 @@ namespace Opm } } + // Initialize saturation scaling parameter + template + void SaturationPropsFromDeck::initEPSKey(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const std::string& keyword, + std::vector& scaleparam) + { + const bool useAqua = phase_usage_.phase_used[Aqua]; + const bool useLiquid = phase_usage_.phase_used[Liquid]; + const bool useVapour = phase_usage_.phase_used[Vapour]; + bool useKeyword = newParserDeck->hasKeyword(keyword); + bool hasENPTVD = newParserDeck->hasKeyword("ENPTVD"); + bool hasENKRVD = newParserDeck->hasKeyword("ENKRVD"); + int itab = 0; + std::vector > > table_dummy; + std::vector > >& table = table_dummy; + + // Active keyword assigned default values for each cell (in case of possible box-wise assignment) + int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + int phase_pos_vapour = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + if ((keyword[0] == 'S' && (useKeyword || hasENPTVD)) || (keyword[1] == 'S' && useKeyword) ) { + if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) { + if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 0))) { + itab = 1; + scaleparam.resize(grid.number_of_cells); + for (int i=0; i 0) { + Opm::EnptvdTable enptvd(newParserDeck->getKeyword("ENPTVD")); + table.resize(1); // only one region supported so far + for (unsigned i = 0; i < table.size(); ++i) { + table[i].resize(9); + + for (int k = 0; k < enptvd.numRows(); ++k) { + table[i][0] = enptvd.getDepthColumn(); + table[i][1] = enptvd.getSwcoColumn(); + table[i][2] = enptvd.getSwcritColumn(); + table[i][3] = enptvd.getSwmaxColumn(); + table[i][4] = enptvd.getSgcoColumn(); + table[i][5] = enptvd.getSgcritColumn(); + table[i][6] = enptvd.getSgmaxColumn(); + table[i][7] = enptvd.getSowcritColumn(); + table[i][8] = enptvd.getSogcritColumn(); + } + } + } + } else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) { + if (keyword == std::string("KRW") || keyword == std::string("IKRW") ) { + if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 0))) { + itab = 1; + scaleparam.resize(grid.number_of_cells); + for (int i=0; i 0) { + Opm::EnkrvdTable enkrvd(newParserDeck->getKeyword("ENKRVD")); + table.resize(1); // only one region supported so far + for (unsigned i = 0; i < table.size(); ++i) { + table[i].resize(8); + + for (int k = 0; k < enkrvd.numRows(); ++k) { + table[i][0] = enkrvd.getDepthColumn(); + table[i][1] = enkrvd.getKrwmaxColumn(); + table[i][2] = enkrvd.getKrgmaxColumn(); + table[i][3] = enkrvd.getKromaxColumn(); + table[i][4] = enkrvd.getKrwcritColumn(); + table[i][5] = enkrvd.getKrgcritColumn(); + table[i][6] = enkrvd.getKrocritgColumn(); + table[i][7] = enkrvd.getKrocritwColumn(); + } + } + } + } + + if (scaleparam.empty()) { + return; + } else if (useKeyword) { + // Keyword values from deck + std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl; + const int* gc = grid.global_cell; + const std::vector& val = newParserDeck->getKeyword(keyword)->getSIDoubleData(); + for (int c = 0; c < int(scaleparam.size()); ++c) { + const int deck_pos = (gc == NULL) ? c : gc[c]; + scaleparam[c] = val[deck_pos]; + } + } else { + // TODO for new parser + /* + std::cout << "--- Scaling parameter '" << keyword << "' assigned via "; + if (keyword[0] == 'S') + newParserDeck.getENPTVD().write(std::cout); + else + newParserDeck.getENKRVD().write(std::cout); + */ + const double* cc = grid.cell_centroids; + const int dim = grid.dimensions; + for (int cell = 0; cell < grid.number_of_cells; ++cell) { + 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]; + if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval + scaleparam[cell] = linearInterpolation(depth, val, zc); + } + } + } + } + } + // Saturation scaling template void SaturationPropsFromDeck::initEPSParam(const int cell, @@ -646,6 +1103,13 @@ namespace Opm { if (scr.empty() && su.empty() && (sxcr.empty() || !do_3pt_) && s0.empty()) { data.doNotScale = true; + data.smin = sl_tab; + if (oil) { + data.smax = (s0_tab < 0.0) ? 1.0 - su_tab : 1.0 - su_tab - s0_tab; + } else { + data.smax = su_tab; + } + data.scr = scr_tab; } else { data.doNotScale = false; data.do_3pt = do_3pt_; diff --git a/opm/core/props/satfunc/SaturationPropsInterface.hpp b/opm/core/props/satfunc/SaturationPropsInterface.hpp index ad2a9439f..04f6730e3 100644 --- a/opm/core/props/satfunc/SaturationPropsInterface.hpp +++ b/opm/core/props/satfunc/SaturationPropsInterface.hpp @@ -73,6 +73,14 @@ namespace Opm const int* cells, double* smin, double* smax) const = 0; + + /// Update saturation state for the hysteresis tracking + /// \param[in] n Number of data points. + /// \param[in] s Array of nP saturation values. + virtual void updateSatHyst(const int n, + const int* cells, + const double* s) = 0; + }; diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 51d2c3b98..090676ee5 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -23,6 +23,7 @@ #include #include #include +#include namespace Opm { @@ -45,28 +46,55 @@ namespace Opm bhp_.resize(nw); wellrates_.resize(nw * np, 0.0); for (int w = 0; w < nw; ++w) { + assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER)); const WellControls* ctrl = wells->ctrls[w]; - // Initialize bhp to be target pressure if - // bhp-controlled well, otherwise set to a little - // above or below (depending on if the well is an - // injector or producer) pressure in first perforation - // cell. - if ((well_controls_get_current(ctrl) < 0) || // SHUT - (well_controls_get_current_type(ctrl) != BHP)) { - const int first_cell = wells->well_cells[wells->well_connpos[w]]; - const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99; - bhp_[w] = safety_factor*state.pressure()[first_cell]; - } else { - bhp_[w] = well_controls_get_current_target( ctrl ); - } - - // Initialize well rates to match controls if type is SURFACE_RATE - if ((well_controls_get_current(ctrl) >= 0) && // open well - (well_controls_get_current_type(ctrl) == SURFACE_RATE)) { - const double rate_target = well_controls_get_current_target(ctrl); - const double * distr = well_controls_get_current_distr( ctrl ); + if (well_controls_well_is_shut(ctrl)) { + // Shut well: + // 1. Assign zero well rates. for (int p = 0; p < np; ++p) { - wellrates_[np*w + p] = rate_target * distr[p]; + wellrates_[np*w + p] = 0.0; + } + // 2. Assign bhp equal to bhp control, if + // applicable, otherwise assign equal to + // first perforation cell pressure. + if (well_controls_get_current_type(ctrl) == BHP) { + bhp_[w] = well_controls_get_current_target( ctrl ); + } else { + const int first_cell = wells->well_cells[wells->well_connpos[w]]; + bhp_[w] = state.pressure()[first_cell]; + } + } else { + // Open well: + // 1. Initialize well rates to match controls + // if type is SURFACE_RATE. Otherwise, we + // cannot set the correct value here, so we + // assign a small rate with the correct + // sign so that any logic depending on that + // sign will work as expected. + if (well_controls_get_current_type(ctrl) == SURFACE_RATE) { + const double rate_target = well_controls_get_current_target(ctrl); + const double * distr = well_controls_get_current_distr( ctrl ); + for (int p = 0; p < np; ++p) { + wellrates_[np*w + p] = rate_target * distr[p]; + } + } else { + const double small_rate = 1e-14; + const double sign = (wells->type[w] == INJECTOR) ? 1.0 : -1.0; + for (int p = 0; p < np; ++p) { + wellrates_[np*w + p] = small_rate * sign; + } + } + // 2. Initialize bhp to be target pressure if + // bhp-controlled well, otherwise set to a + // little above or below (depending on if + // the well is an injector or producer) + // pressure in first perforation cell. + if (well_controls_get_current_type(ctrl) == BHP) { + bhp_[w] = well_controls_get_current_target( ctrl ); + } else { + const int first_cell = wells->well_cells[wells->well_connpos[w]]; + const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99; + bhp_[w] = safety_factor*state.pressure()[first_cell]; } } } diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index 12dbe7b52..89ff89234 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -32,6 +32,8 @@ #include #include +#include + #include #include @@ -471,6 +473,108 @@ namespace Opm } + /// Initialize a state from input deck. + template + void initStateFromDeck(const UnstructuredGrid& grid, + const Props& props, + Opm::DeckConstPtr newParserDeck, + const double gravity, + State& state) + { + const int num_phases = props.numPhases(); + const PhaseUsage pu = phaseUsageFromDeck(newParserDeck); + if (num_phases != pu.num_phases) { + 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); + if (newParserDeck->hasKeyword("EQUIL") && newParserDeck->hasKeyword("PRESSURE")) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial " + "condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)"); + } + + if (newParserDeck->hasKeyword("EQUIL")) { + if (num_phases != 2) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios."); + } + if (pu.phase_used[BlackoilPhases::Vapour]) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas)."); + } + // Set saturations depending on oil-water contact. + EquilWrapper equil(newParserDeck->getKeyword("EQUIL")); + if (equil.numRegions() != 1) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet."); + } + const double woc = equil.waterOilContactDepth(0); + initWaterOilContact(grid, props, woc, WaterBelow, state); + // Set pressure depending on densities and depths. + const double datum_z = equil.datumDepth(0); + const double datum_p = equil.datumDepthPressure(0); + initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state); + } else if (newParserDeck->hasKeyword("PRESSURE")) { + // Set saturations from SWAT/SGAS, pressure from PRESSURE. + std::vector& s = state.saturation(); + std::vector& p = state.pressure(); + const std::vector& p_deck = newParserDeck->getKeyword("PRESSURE")->getSIDoubleData(); + const int num_cells = grid.number_of_cells; + if (num_phases == 2) { + if (!pu.phase_used[BlackoilPhases::Aqua]) { + // oil-gas: we require SGAS + if (!newParserDeck->hasKeyword("SGAS")) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS keyword in 2-phase init"); + } + const std::vector& sg_deck = newParserDeck->getKeyword("SGAS")->getSIDoubleData(); + 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]; + s[2*c + gpos] = sg_deck[c_deck]; + s[2*c + opos] = 1.0 - sg_deck[c_deck]; + p[c] = p_deck[c_deck]; + } + } else { + // water-oil or water-gas: we require SWAT + if (!newParserDeck->hasKeyword("SWAT")) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SWAT keyword in 2-phase init"); + } + const std::vector& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData(); + 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]; + s[2*c + wpos] = sw_deck[c_deck]; + s[2*c + nwpos] = 1.0 - sw_deck[c_deck]; + p[c] = p_deck[c_deck]; + } + } + } else if (num_phases == 3) { + const bool has_swat_sgas = newParserDeck->hasKeyword("SWAT") && newParserDeck->hasKeyword("SGAS"); + if (!has_swat_sgas) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init."); + } + const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const int opos = pu.phase_pos[BlackoilPhases::Liquid]; + const std::vector& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData(); + const std::vector& sg_deck = newParserDeck->getKeyword("SGAS")->getSIDoubleData(); + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.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]; + p[c] = p_deck[c_deck]; + } + } else { + OPM_THROW(std::runtime_error, "initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases."); + } + } else { + OPM_THROW(std::runtime_error, "initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS."); + } + + // Finally, init face pressures. + initFacePressure(grid, state); + } + /// Initialize a state from input deck. template void initStateFromDeck(const UnstructuredGrid& grid, @@ -485,6 +589,10 @@ namespace Opm OPM_THROW(std::runtime_error, "initStateFromDeck(): user specified property object with " << num_phases << " phases, " "found " << pu.num_phases << " phases in deck."); } + if (deck.hasField("EQUIL") && deck.hasField("PRESSURE")) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial " + "condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)"); + } state.init(grid, num_phases); if (deck.hasField("EQUIL")) { if (num_phases != 2) { @@ -568,10 +676,6 @@ namespace Opm initFacePressure(grid, state); } - - - - /// Initialize surface volume from pressure and saturation by z = As. /// Here saturation is used as an intial guess for z in the /// computation of A. @@ -770,6 +874,39 @@ namespace Opm } } + /// Initialize a blackoil state from input deck. + template + void initBlackoilStateFromDeck(const UnstructuredGrid& grid, + const Props& props, + Opm::DeckConstPtr newParserDeck, + const double gravity, + State& state) + { + initStateFromDeck(grid, props, newParserDeck, gravity, state); + if (newParserDeck->hasKeyword("RS")) { + const std::vector& rs_deck = newParserDeck->getKeyword("RS")->getSIDoubleData(); + 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]; + state.gasoilratio()[c] = rs_deck[c_deck]; + } + initBlackoilSurfvolUsingRSorRV(grid, props, state); + computeSaturation(props,state); + } else if (newParserDeck->hasKeyword("RV")){ + const std::vector& rv_deck = newParserDeck->getKeyword("RV")->getSIDoubleData(); + 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]; + state.rv()[c] = rv_deck[c_deck]; + } + initBlackoilSurfvolUsingRSorRV(grid, props, state); + computeSaturation(props,state); + } + else { + OPM_THROW(std::runtime_error, "Temporarily, we require the RS or the RV field."); + } + } + } // namespace Opm diff --git a/opm/core/well_controls.h b/opm/core/well_controls.h index 70ea1b49b..6998f4829 100644 --- a/opm/core/well_controls.h +++ b/opm/core/well_controls.h @@ -34,10 +34,10 @@ enum WellControlType { SURFACE_RATE /**< Well constrained by surface volume flow rate */ }; -struct WellControls; +struct WellControls; bool -well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2); +well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2 , bool verbose); struct WellControls * well_controls_create(void); @@ -55,8 +55,18 @@ well_controls_get_current( const struct WellControls * ctrl); void well_controls_set_current( struct WellControls * ctrl, int current); -void -well_controls_invert_current( struct WellControls * ctrl ); + +bool +well_controls_well_is_shut(const struct WellControls * ctrl); + +bool +well_controls_well_is_open(const struct WellControls * ctrl); + +void +well_controls_open_well( struct WellControls * ctrl); + +void +well_controls_shut_well( struct WellControls * ctrl); int well_controls_add_new(enum WellControlType type , double target , const double * distr , struct WellControls * ctrl); diff --git a/opm/core/wells.h b/opm/core/wells.h index dfc9be176..05b6ce3e8 100644 --- a/opm/core/wells.h +++ b/opm/core/wells.h @@ -263,7 +263,7 @@ struct Wells * clone_wells(const struct Wells *W); bool -wells_equal(const struct Wells *W1, const struct Wells *W2); +wells_equal(const struct Wells *W1, const struct Wells *W2 , bool verbose); diff --git a/opm/core/wells/WellCollection.cpp b/opm/core/wells/WellCollection.cpp index 350d2ffb3..9218db06e 100644 --- a/opm/core/wells/WellCollection.cpp +++ b/opm/core/wells/WellCollection.cpp @@ -19,10 +19,66 @@ #include "config.h" #include + +#include +#include + +#include + #include namespace Opm { + void WellCollection::addField(GroupConstPtr fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage) { + WellsGroupInterface* fieldNode = findNode(fieldGroup->name()); + if (fieldNode) { + OPM_THROW(std::runtime_error, "Trying to add FIELD node, but this already exists. Can only have one FIELD node."); + } + + roots_.push_back(createGroupWellsGroup(fieldGroup, timeStep, phaseUsage)); + } + + void WellCollection::addGroup(GroupConstPtr groupChild, std::string parent_name, + size_t timeStep, const PhaseUsage& phaseUsage) { + WellsGroupInterface* parent = findNode(parent_name); + if (!parent) { + OPM_THROW(std::runtime_error, "Trying to add child group to group named " << parent_name << ", but this does not exist in the WellCollection."); + } + + if (findNode(groupChild->name())) { + OPM_THROW(std::runtime_error, "Trying to add child group named " << groupChild->name() << ", but this group is already in the WellCollection."); + + } + + std::shared_ptr child = createGroupWellsGroup(groupChild, timeStep, phaseUsage); + + WellsGroup* parent_as_group = static_cast (parent); + if (!parent_as_group) { + OPM_THROW(std::runtime_error, "Trying to add child group to group named " << parent->name() << ", but it's not a group."); + } + parent_as_group->addChild(child); + child->setParent(parent); + } + + void WellCollection::addWell(WellConstPtr wellChild, size_t timeStep, const PhaseUsage& phaseUsage) { + WellsGroupInterface* parent = findNode(wellChild->getGroupName(timeStep)); + if (!parent) { + OPM_THROW(std::runtime_error, "Trying to add well " << wellChild->name() << " Step: " << boost::lexical_cast(timeStep) << " to group named " << wellChild->getGroupName(timeStep) << ", but this group does not exist in the WellCollection."); + } + + std::shared_ptr child = createWellWellsGroup(wellChild, timeStep, phaseUsage); + + WellsGroup* parent_as_group = static_cast (parent); + if (!parent_as_group) { + OPM_THROW(std::runtime_error, "Trying to add well to group named " << wellChild->getGroupName(timeStep) << ", but it's not a group."); + } + parent_as_group->addChild(child); + + leaf_nodes_.push_back(static_cast(child.get())); + + child->setParent(parent); + } + void WellCollection::addChild(const std::string& child_name, const std::string& parent_name, @@ -33,6 +89,7 @@ namespace Opm roots_.push_back(createWellsGroup(parent_name, deck)); parent = roots_[roots_.size() - 1].get(); } + std::shared_ptr child; for (size_t i = 0; i < roots_.size(); ++i) { @@ -47,6 +104,7 @@ namespace Opm break; } } + if (!child.get()) { child = createWellsGroup(child_name, deck); } @@ -65,7 +123,6 @@ namespace Opm } - const std::vector& WellCollection::getLeafNodes() const { return leaf_nodes_; } diff --git a/opm/core/wells/WellCollection.hpp b/opm/core/wells/WellCollection.hpp index 9cdb9e7f1..f084a04d8 100644 --- a/opm/core/wells/WellCollection.hpp +++ b/opm/core/wells/WellCollection.hpp @@ -23,10 +23,15 @@ #define OPM_WELLCOLLECTION_HPP #include +#include + #include #include #include -#include +#include + +#include +#include namespace Opm { @@ -34,6 +39,14 @@ namespace Opm class WellCollection { public: + + void addField(GroupConstPtr fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage); + + void addWell(WellConstPtr wellChild, size_t timeStep, const PhaseUsage& phaseUsage); + + void addGroup(GroupConstPtr groupChild, std::string parent_name, + size_t timeStep, const PhaseUsage& phaseUsage); + /// Adds and creates if necessary the child to the collection /// and appends it to parent's children. Also adds and creates the parent /// if necessary. diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index 6c1059a29..050e18712 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -74,7 +74,7 @@ namespace Opm { return parent_; } - const std::string& WellsGroupInterface::name() + const std::string& WellsGroupInterface::name() const { return name_; } @@ -747,9 +747,7 @@ namespace Opm void WellNode::shutWell() { if (shut_well_) { - // We set the tilde of the current control - // set_current_control(self_index_, -1, wells_); - well_controls_invert_current(wells_->ctrls[self_index_]); + well_controls_shut_well( wells_->ctrls[self_index_]); } else { const double target = 0.0; @@ -767,7 +765,7 @@ namespace Opm well_controls_iset_target( wells_->ctrls[self_index_] , group_control_index_ , target); well_controls_iset_distr(wells_->ctrls[self_index_] , group_control_index_ , distr); } - well_controls_invert_current(wells_->ctrls[self_index_]); + well_controls_open_well( wells_->ctrls[self_index_]); } } @@ -1141,4 +1139,56 @@ namespace Opm return return_value; } + + std::shared_ptr createGroupWellsGroup(GroupConstPtr group, size_t timeStep, const PhaseUsage& phase_usage ) + { + InjectionSpecification injection_specification; + ProductionSpecification production_specification; + if (group->isInjectionGroup(timeStep)) { + injection_specification.injector_type_ = toInjectorType(Phase::PhaseEnum2String(group->getInjectionPhase(timeStep))); + injection_specification.control_mode_ = toInjectionControlMode(GroupInjection::ControlEnum2String(group->getInjectionControlMode(timeStep))); + injection_specification.surface_flow_max_rate_ = group->getSurfaceMaxRate(timeStep); + injection_specification.reservoir_flow_max_rate_ = group->getReservoirMaxRate(timeStep); + injection_specification.reinjection_fraction_target_ = group->getTargetReinjectFraction(timeStep); + injection_specification.voidage_replacment_fraction_ = group->getTargetVoidReplacementFraction(timeStep); + } + else if (group->isProductionGroup(timeStep)) { + production_specification.oil_max_rate_ = group->getOilTargetRate(timeStep); + production_specification.control_mode_ = toProductionControlMode(GroupProduction::ControlEnum2String(group->getProductionControlMode(timeStep))); + production_specification.water_max_rate_ = group->getWaterTargetRate(timeStep); + production_specification.gas_max_rate_ = group->getGasTargetRate(timeStep); + production_specification.liquid_max_rate_ = group->getLiquidTargetRate(timeStep); + production_specification.procedure_ = toProductionProcedure(GroupProductionExceedLimit::ActionEnum2String(group->getProductionExceedLimitAction(timeStep))); + production_specification.reservoir_flow_max_rate_ = group->getReservoirMaxRate(timeStep); + } + + std::shared_ptr wells_group(new WellsGroup(group->name(), production_specification, injection_specification, phase_usage)); + return wells_group; + } + + std::shared_ptr createWellWellsGroup(WellConstPtr well, size_t timeStep, const PhaseUsage& phase_usage ) + { + InjectionSpecification injection_specification; + ProductionSpecification production_specification; + if (well->isInjector(timeStep)) { + const WellInjectionProperties& properties = well->getInjectionProperties(timeStep); + injection_specification.BHP_limit_ = properties.BHPLimit; + injection_specification.injector_type_ = toInjectorType(WellInjector::Type2String(properties.injectorType)); + injection_specification.control_mode_ = toInjectionControlMode(WellInjector::ControlMode2String(properties.controlMode)); + injection_specification.surface_flow_max_rate_ = properties.surfaceInjectionRate; + injection_specification.reservoir_flow_max_rate_ = properties.reservoirInjectionRate; + production_specification.guide_rate_ = 0.0; // We know we're not a producer + } + else if (well->isProducer(timeStep)) { + const WellProductionProperties& properties = well->getProductionProperties(timeStep); + production_specification.BHP_limit_ = properties.BHPLimit; + production_specification.reservoir_flow_max_rate_ = properties.ResVRate; + production_specification.oil_max_rate_ = properties.OilRate; + production_specification.control_mode_ = toProductionControlMode(WellProducer::ControlMode2String(properties.controlMode)); + production_specification.water_max_rate_ = properties.WaterRate; + injection_specification.guide_rate_ = 0.0; // we know we're not an injector + } + std::shared_ptr wells_group(new WellNode(well->name(), production_specification, injection_specification, phase_usage)); + return wells_group; + } } diff --git a/opm/core/wells/WellsGroup.hpp b/opm/core/wells/WellsGroup.hpp index 143da9357..03a61d5e1 100644 --- a/opm/core/wells/WellsGroup.hpp +++ b/opm/core/wells/WellsGroup.hpp @@ -25,6 +25,10 @@ #include #include #include + +#include +#include + #include #include @@ -58,7 +62,7 @@ namespace Opm virtual ~WellsGroupInterface(); /// The unique identifier for the well or well group. - const std::string& name(); + const std::string& name() const; /// Production specifications for the well or well group. const ProductionSpecification& prodSpec() const; @@ -403,9 +407,21 @@ namespace Opm /// \param[in] name the name of the wells group. /// \param[in] deck the deck from which to fetch information. std::shared_ptr createWellsGroup(const std::string& name, - const EclipseGridParser& deck); + const EclipseGridParser& deck); + /// Creates the WellsGroupInterface for the given well + /// \param[in] well the Well to construct object for + /// \param[in] timeStep the time step in question + /// \param[in] the phase usage + std::shared_ptr createWellWellsGroup(WellConstPtr well, size_t timeStep, + const PhaseUsage& phase_usage ); + /// Creates the WellsGroupInterface for the given Group + /// \param[in] group the Group to construct object for + /// \param[in] timeStep the time step in question + /// \param[in] the phase usage + std::shared_ptr createGroupWellsGroup(GroupConstPtr group, size_t timeStep, + const PhaseUsage& phase_usage ); } #endif /* OPM_WELLSGROUP_HPP */ diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 007c8b52c..02e762094 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -18,6 +18,8 @@ */ #include "config.h" + + #include #include #include @@ -28,6 +30,8 @@ #include #include +#include + #include #include #include @@ -43,23 +47,8 @@ namespace { - struct WellData - { - WellType type; - // WellControlType control; - // double target; - double reference_bhp_depth; - // Opm::InjectionSpecification::InjectorType injected_phase; - int welspecsline; - }; - struct PerfData - { - int cell; - double well_index; - }; - namespace ProductionControl { enum Mode { ORAT, WRAT, GRAT, @@ -101,6 +90,33 @@ namespace << control << " in input file"); } } + + + Mode mode(Opm::WellProducer::ControlModeEnum controlMode) + { + switch( controlMode ) { + case Opm::WellProducer::ORAT: + return ORAT; + case Opm::WellProducer::WRAT: + return WRAT; + case Opm::WellProducer::GRAT: + return GRAT; + case Opm::WellProducer::LRAT: + return LRAT; + case Opm::WellProducer::CRAT: + return CRAT; + case Opm::WellProducer::RESV: + return RESV; + case Opm::WellProducer::BHP: + return BHP; + case Opm::WellProducer::THP: + return THP; + case Opm::WellProducer::GRUP: + return GRUP; + default: + throw std::invalid_argument("unhandled enum value"); + } + } } // namespace ProductionControl @@ -140,6 +156,25 @@ namespace << control << " in input file"); } } + + Mode mode(Opm::WellInjector::ControlModeEnum controlMode) + { + switch ( controlMode ) { + case Opm::WellInjector::GRUP: + return GRUP; + case Opm::WellInjector::RESV: + return RESV; + case Opm::WellInjector::RATE: + return RATE; + case Opm::WellInjector::THP: + return THP; + case Opm::WellInjector::BHP: + return BHP; + default: + throw std::invalid_argument("unhandled enum value"); + } + } + } // namespace InjectionControl @@ -233,10 +268,9 @@ namespace Opm { } - - /// Construct wells from deck. - WellsManager::WellsManager(const Opm::EclipseGridParser& deck, + WellsManager::WellsManager(const Opm::EclipseStateConstPtr eclipseState, + const size_t timeStep, const UnstructuredGrid& grid, const double* permeability) : w_(0) @@ -244,484 +278,52 @@ namespace Opm if (grid.dimensions != 3) { OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional."); } - // NOTE: Implementation copied and modified from dune-porsol's class BlackoilWells. - std::vector keywords; - keywords.push_back("WELSPECS"); - keywords.push_back("COMPDAT"); -// keywords.push_back("WELTARG"); - if (!deck.hasFields(keywords)) { - OPM_MESSAGE("Missing well keywords in deck, initializing no wells."); + + if (eclipseState->getSchedule()->numWells() == 0) { + OPM_MESSAGE("No wells specified in Schedule section, initializing no wells"); return; } - if (!(deck.hasField("WCONINJE") || deck.hasField("WCONPROD")) ) { - OPM_THROW(std::runtime_error, "Needed field is missing in file"); - } + + std::map cartesian_to_compressed; + setupCompressedToCartesian(grid, cartesian_to_compressed); // Obtain phase usage data. - PhaseUsage pu = phaseUsageFromDeck(deck); + PhaseUsage pu = phaseUsageFromDeck(eclipseState); // 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); + ScheduleConstPtr schedule = eclipseState->getSchedule(); + std::vector wells = schedule->getWells(timeStep); + well_names.reserve(wells.size()); + well_data.reserve(wells.size()); - // global_cell is a map from compressed cells to Cartesian grid cells. - // We must make the inverse lookup. - const int* global_cell = grid.global_cell; - const int* cpgdim = grid.cartdims; - std::map cartesian_to_compressed; + createWellsFromSpecs(wells, timeStep, grid, well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability); + setupWellControls(wells, timeStep, well_names, pu); - if (global_cell) { - for (int i = 0; i < grid.number_of_cells; ++i) { - cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); - } - } - else { - for (int i = 0; i < grid.number_of_cells; ++i) { - cartesian_to_compressed.insert(std::make_pair(i, i)); - } + { + GroupTreeNodeConstPtr fieldNode = eclipseState->getSchedule()->getGroupTree(timeStep)->getNode("FIELD"); + GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(fieldNode->name()); + well_collection_.addField(fieldGroup, timeStep, pu); + addChildGroups(fieldNode, eclipseState->getSchedule(), timeStep, pu); } - // Get COMPDAT data - // It is *not* allowed to have multiple lines corresponding to - // the same perforation! - const COMPDAT& compdat = deck.getCOMPDAT(); - const int num_compdat = compdat.compdat.size(); - for (int kw = 0; kw < num_compdat; ++kw) { - // Extract well name, or the part of the well name that - // comes before the '*'. - std::string name = compdat.compdat[kw].well_; - std::string::size_type len = name.find('*'); - if (len != std::string::npos) { - name = name.substr(0, len); - } - // Look for well with matching name. - bool found = false; - for (int wix = 0; wix < num_wells; ++wix) { - if (well_names[wix].compare(0,len, name) == 0) { // equal - // Extract corresponding WELSPECS defintion for - // purpose of default location specification. - const WelspecsLine& wspec = welspecs.welspecs[well_data[wix].welspecsline]; - - // We have a matching name. - int ix = compdat.compdat[kw].grid_ind_[0] - 1; - int jy = compdat.compdat[kw].grid_ind_[1] - 1; - int kz1 = compdat.compdat[kw].grid_ind_[2] - 1; - int kz2 = compdat.compdat[kw].grid_ind_[3] - 1; - - if (ix < 0) { - // Defaulted I location. Extract from WELSPECS. - ix = wspec.I_ - 1; - } - if (jy < 0) { - // Defaulted J location. Extract from WELSPECS. - jy = wspec.J_ - 1; - } - if (kz1 < 0) { - // Defaulted KZ1. Use top layer. - kz1 = 0; - } - if (kz2 < 0) { - // Defaulted KZ2. Use bottom layer. - kz2 = cpgdim[2] - 1; - } - - for (int kz = kz1; kz <= kz2; ++kz) { - int cart_grid_indx = ix + cpgdim[0]*(jy + cpgdim[1]*kz); - std::map::const_iterator cgit = - cartesian_to_compressed.find(cart_grid_indx); - if (cgit == cartesian_to_compressed.end()) { - OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << ix << ' ' << jy << ' ' - << kz << " not found in grid (well = " << name << ')'); - } - int cell = cgit->second; - PerfData pd; - pd.cell = cell; - if (compdat.compdat[kw].connect_trans_fac_ > 0.0) { - pd.well_index = compdat.compdat[kw].connect_trans_fac_; - } else { - double radius = 0.5*compdat.compdat[kw].diameter_; - if (radius <= 0.0) { - radius = 0.5*unit::feet; - OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); - } - std::array cubical = getCubeDim(grid, cell); - const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell]; - pd.well_index = computeWellIndex(radius, cubical, cell_perm, - compdat.compdat[kw].skin_factor_); - } - wellperf_data[wix].push_back(pd); - } - found = true; - break; - } - } - if (!found) { - OPM_THROW(std::runtime_error, "Undefined well name: " << compdat.compdat[kw].well_ - << " in COMPDAT"); - } + for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) { + well_collection_.addWell((*wellIter), timeStep, pu); } - // Set up reference depths that were defaulted. Count perfs. - int num_perfs = 0; - assert(grid.dimensions == 3); - for (int w = 0; w < num_wells; ++w) { - num_perfs += wellperf_data[w].size(); - if (well_data[w].reference_bhp_depth < 0.0) { - // It was defaulted. Set reference depth to minimum perforation depth. - double min_depth = 1e100; - int num_wperfs = wellperf_data[w].size(); - for (int perf = 0; perf < num_wperfs; ++perf) { - double depth = grid.cell_centroids[3*wellperf_data[w][perf].cell + 2]; - min_depth = std::min(min_depth, depth); - } - well_data[w].reference_bhp_depth = min_depth; - } - } + well_collection_.setWellsPointer(w_); + well_collection_.applyGroupControls(); - // 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"); - } - } - */ - } + setupGuideRates(wells, timeStep, well_data, well_names_to_index); // Debug output. #define EXTRA_OUTPUT @@ -743,89 +345,13 @@ namespace Opm } */ #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(); } + + + /// Destructor. WellsManager::~WellsManager() { @@ -879,4 +405,408 @@ namespace Opm well_collection_.applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase); } + void WellsManager::setupCompressedToCartesian(const UnstructuredGrid& grid, 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) { + cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); + } + } + else { + for (int i = 0; i < grid.number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(i, i)); + } + } + + } + + void WellsManager::createWellsFromSpecs(std::vector& wells, size_t timeStep, + const UnstructuredGrid& grid, + std::vector& well_names, + std::vector& well_data, + std::map& well_names_to_index, + const PhaseUsage& phaseUsage, + std::map cartesian_to_compressed, + const double* permeability) + { + std::vector > wellperf_data; + wellperf_data.resize(wells.size()); + + int well_index = 0; + for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) { + WellConstPtr well = (*wellIter); + { // WELSPECS handling + well_names_to_index[well->name()] = well_index; + well_names.push_back(well->name()); + { + WellData wd; + // If negative (defaulted), set refdepth to a marker + // value, will be changed after getting perforation + // data to the centroid of the cell of the top well + // perforation. + wd.reference_bhp_depth = (well->getRefDepth() < 0.0) ? -1e100 : well->getRefDepth(); + wd.welspecsline = -1; + if (well->isInjector( timeStep )) + wd.type = INJECTOR; + else + wd.type = PRODUCER; + well_data.push_back(wd); + } + } + + { // COMPDAT handling + CompletionSetConstPtr completionSet = well->getCompletions(timeStep); + for (size_t c=0; csize(); c++) { + CompletionConstPtr completion = completionSet->get(c); + int i = completion->getI(); + int j = completion->getJ(); + int k = completion->getK(); + + const int* cpgdim = grid.cartdims; + int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k); + std::map::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx); + if (cgit == cartesian_to_compressed.end()) { + OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' ' + << k << " not found in grid (well = " << well->name() << ')'); + } + int cell = cgit->second; + PerfData pd; + pd.cell = cell; + if (completion->getCF() > 0.0) { + pd.well_index = completion->getCF(); + } else { + double radius = 0.5*completion->getDiameter(); + if (radius <= 0.0) { + radius = 0.5*unit::feet; + OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); + } + std::array cubical = getCubeDim(grid, cell); + const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell]; + pd.well_index = computeWellIndex(radius, cubical, cell_perm, completion->getDiameter()); + } + wellperf_data[well_index].push_back(pd); + } + } + well_index++; + } + + // Set up reference depths that were defaulted. Count perfs. + + const int num_wells = well_data.size(); + + int num_perfs = 0; + assert(grid.dimensions == 3); + for (int w = 0; w < num_wells; ++w) { + num_perfs += wellperf_data[w].size(); + if (well_data[w].reference_bhp_depth < 0.0) { + // It was defaulted. Set reference depth to minimum perforation depth. + double min_depth = 1e100; + int num_wperfs = wellperf_data[w].size(); + for (int perf = 0; perf < num_wperfs; ++perf) { + double depth = grid.cell_centroids[3*wellperf_data[w][perf].cell + 2]; + min_depth = std::min(min_depth, depth); + } + well_data[w].reference_bhp_depth = min_depth; + } + } + + // Create the well data structures. + w_ = create_wells(phaseUsage.num_phases, num_wells, num_perfs); + if (!w_) { + OPM_THROW(std::runtime_error, "Failed creating Wells struct."); + } + + + // Add wells. + for (int w = 0; w < num_wells; ++w) { + const int w_num_perf = wellperf_data[w].size(); + std::vector perf_cells(w_num_perf); + std::vector perf_prodind(w_num_perf); + for (int perf = 0; perf < w_num_perf; ++perf) { + perf_cells[perf] = wellperf_data[w][perf].cell; + perf_prodind[perf] = wellperf_data[w][perf].well_index; + } + const double* comp_frac = NULL; + // We initialize all wells with a null component fraction, + // and must (for injection wells) overwrite it later. + int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, w_num_perf, + comp_frac, &perf_cells[0], &perf_prodind[0], well_names[w].c_str(), w_); + if (!ok) { + OPM_THROW(std::runtime_error, "Failed adding well " << well_names[w] << " to Wells data structure."); + } + } + + } + + void WellsManager::setupWellControls(std::vector& wells, size_t timeStep, + std::vector& well_names, const PhaseUsage& phaseUsage) { + int well_index = 0; + for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) { + WellConstPtr well = (*wellIter); + + if ( !( well->getStatus( timeStep ) == WellCommon::SHUT || well->getStatus( timeStep ) == WellCommon::OPEN) ) { + OPM_THROW(std::runtime_error, "Currently we do not support well status " << WellCommon::Status2String(well->getStatus( timeStep ))); + } + + if (well->isInjector(timeStep)) { + const WellInjectionProperties& injectionProperties = well->getInjectionProperties(timeStep); + int ok = 1; + int control_pos[5] = { -1, -1, -1, -1, -1 }; + + clear_well_controls(well_index, w_); + if (injectionProperties.hasInjectionControl(WellInjector::RATE)) { + control_pos[InjectionControl::RATE] = well_controls_get_num(w_->ctrls[well_index]); + double distr[3] = { 0.0, 0.0, 0.0 }; + WellInjector::TypeEnum injectorType = injectionProperties.injectorType; + + if (injectorType == WellInjector::TypeEnum::WATER) { + distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0; + } else if (injectorType == WellInjector::TypeEnum::OIL) { + distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0; + } else if (injectorType == WellInjector::TypeEnum::GAS) { + distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0; + } + + ok = append_well_controls(SURFACE_RATE, + injectionProperties.surfaceInjectionRate, + distr, + well_index, + w_); + } + + if (ok && injectionProperties.hasInjectionControl(WellInjector::RESV)) { + control_pos[InjectionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]); + double distr[3] = { 0.0, 0.0, 0.0 }; + WellInjector::TypeEnum injectorType = injectionProperties.injectorType; + + if (injectorType == WellInjector::TypeEnum::WATER) { + distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0; + } else if (injectorType == WellInjector::TypeEnum::OIL) { + distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0; + } else if (injectorType == WellInjector::TypeEnum::GAS) { + distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0; + } + + ok = append_well_controls(RESERVOIR_RATE, + injectionProperties.reservoirInjectionRate, + distr, + well_index, + w_); + } + + if (ok && injectionProperties.hasInjectionControl(WellInjector::BHP)) { + control_pos[InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]); + ok = append_well_controls(BHP, + injectionProperties.BHPLimit, + NULL, + well_index, + w_); + } + + if (ok && injectionProperties.hasInjectionControl(WellInjector::THP)) { + OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[well_index]); + } + + + if (!ok) { + OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[well_index]); + } + + + { + InjectionControl::Mode mode = InjectionControl::mode( injectionProperties.controlMode ); + int cpos = control_pos[mode]; + if (cpos == -1 && mode != InjectionControl::GRUP) { + OPM_THROW(std::runtime_error, "Control not specified in well " << well_names[well_index]); + } + + // We need to check if the well is shut or not + if (well->getStatus( timeStep ) == WellCommon::SHUT) { + well_controls_shut_well( w_->ctrls[well_index] ); + } + set_current_control(well_index, cpos, w_); + } + + + // Set well component fraction. + double cf[3] = { 0.0, 0.0, 0.0 }; + { + WellInjector::TypeEnum injectorType = injectionProperties.injectorType; + + if (injectorType == WellInjector::WATER) { + if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) { + OPM_THROW(std::runtime_error, "Water phase not used, yet found water-injecting well."); + } + cf[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0; + } else if (injectorType == WellInjector::OIL) { + if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) { + OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-injecting well."); + } + cf[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0; + } else if (injectorType == WellInjector::GAS) { + if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) { + OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-injecting well."); + } + cf[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0; + } + std::copy(cf, cf + phaseUsage.num_phases, w_->comp_frac + well_index*phaseUsage.num_phases); + } + } + + if (well->isProducer(timeStep)) { + // Add all controls that are present in well. + // First we must clear existing controls, in case the + // current WCONPROD line is modifying earlier controls. + const WellProductionProperties& productionProperties = well->getProductionProperties(timeStep); + int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 }; + int ok = 1; + + clear_well_controls(well_index, w_); + if (ok && productionProperties.hasProductionControl(WellProducer::ORAT)) { + if (!phaseUsage.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[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, + -productionProperties.OilRate, + distr, + well_index, + w_); + } + + if (ok && productionProperties.hasProductionControl(WellProducer::WRAT)) { + 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]); + double distr[3] = { 0.0, 0.0, 0.0 }; + distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0; + ok = append_well_controls(SURFACE_RATE, + -productionProperties.WaterRate, + distr, + well_index, + w_); + } + + if (ok && productionProperties.hasProductionControl(WellProducer::GRAT)) { + 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]); + double distr[3] = { 0.0, 0.0, 0.0 }; + distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0; + ok = append_well_controls(SURFACE_RATE, + -productionProperties.GasRate, + distr, + well_index, + w_); + } + + if (ok && productionProperties.hasProductionControl(WellProducer::LRAT)) { + if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) { + OPM_THROW(std::runtime_error, "Water phase not active and LRAT control specified."); + } + 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]); + 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; + ok = append_well_controls(SURFACE_RATE, + -productionProperties.LiquidRate , + distr, + well_index, + w_); + } + + if (ok && productionProperties.hasProductionControl(WellProducer::RESV)) { + control_pos[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, + -productionProperties.ResVRate , + distr, + well_index, + w_); + } + + if (ok && productionProperties.hasProductionControl(WellProducer::BHP)) { + control_pos[ProductionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]); + ok = append_well_controls(BHP, + productionProperties.BHPLimit , + NULL, + well_index, + w_); + } + + if (ok && productionProperties.hasProductionControl(WellProducer::THP)) { + OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[well_index]); + } + + if (!ok) { + OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[well_index]); + } + + ProductionControl::Mode mode = ProductionControl::mode(productionProperties.controlMode); + int cpos = control_pos[mode]; + if (well->getStatus(timeStep) == WellCommon::SHUT) { + well_controls_shut_well( w_->ctrls[well_index] ); + } else if (cpos == -1 && mode != ProductionControl::GRUP) { + OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]); + } + set_current_control(well_index, cpos, w_); + } + well_index++; + } + + } + + void WellsManager::addChildGroups(GroupTreeNodeConstPtr parentNode, ScheduleConstPtr schedule, size_t timeStep, const PhaseUsage& phaseUsage) { + for (auto childIter = parentNode->begin(); childIter != parentNode->end(); ++childIter) { + GroupTreeNodeConstPtr childNode = (*childIter).second; + well_collection_.addGroup(schedule->getGroup(childNode->name()), parentNode->name(), timeStep, phaseUsage); + addChildGroups(childNode, schedule, timeStep, phaseUsage); + } + } + + + + + void WellsManager::setupGuideRates(std::vector& wells, const size_t timeStep, std::vector& well_data, std::map& well_names_to_index) + { + for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) { + WellConstPtr well = *wellIter; + const int wix = well_names_to_index[well->name()]; + WellNode& wellnode = *well_collection_.getLeafNodes()[wix]; + + if (well->getGuideRatePhase(timeStep) != GuideRate::UNDEFINED) { + if (well_data[wix].type == PRODUCER) { + wellnode.prodSpec().guide_rate_ = well->getGuideRate(timeStep); + if (well->getGuideRatePhase(timeStep) == GuideRate::OIL) { + wellnode.prodSpec().guide_rate_type_ = ProductionSpecification::OIL; + } else { + OPM_THROW(std::runtime_error, "Guide rate type " << GuideRate::GuideRatePhaseEnum2String(well->getGuideRatePhase(timeStep)) << " specified for producer " + << well->name() << " in WGRUPCON, cannot handle."); + } + } else if (well_data[wix].type == INJECTOR) { + wellnode.injSpec().guide_rate_ = well->getGuideRate(timeStep); + if (well->getGuideRatePhase(timeStep) == GuideRate::RAT) { + wellnode.injSpec().guide_rate_type_ = InjectionSpecification::RAT; + } else { + OPM_THROW(std::runtime_error, "Guide rate type " << GuideRate::GuideRatePhaseEnum2String(well->getGuideRatePhase(timeStep)) << " specified for injector " + << well->name() << " in WGRUPCON, cannot handle."); + } + } else { + OPM_THROW(std::runtime_error, "Unknown well type " << well_data[wix].type << " for well " << well->name()); + } + } + } + } + } // namespace Opm diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp index 35f7f471d..761a32315 100644 --- a/opm/core/wells/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -21,8 +21,11 @@ #define OPM_WELLSMANAGER_HEADER_INCLUDED +#include + #include #include +#include struct Wells; struct UnstructuredGrid; @@ -32,7 +35,22 @@ namespace Opm { class EclipseGridParser; + struct WellData + { + WellType type; + // WellControlType control; + // double target; + double reference_bhp_depth; + // Opm::InjectionSpecification::InjectorType injected_phase; + int welspecsline; + }; + + struct PerfData + { + int cell; + double well_index; + }; /// This class manages a Wells struct in the sense that it /// encapsulates creation and destruction of the wells /// data structure. @@ -40,34 +58,35 @@ namespace Opm class WellsManager { public: - /// Default constructor -- no wells. - WellsManager(); + /// Default constructor -- no wells. + WellsManager(); - /// Construct from existing wells object. - /// WellsManager is not properly initialised in the sense that the logic to - /// manage control switching does not exist. - /// - /// @param[in] W Existing wells object. - WellsManager(struct Wells* W); + /// Construct from existing wells object. + /// WellsManager is not properly initialised in the sense that the logic to + /// manage control switching does not exist. + /// + /// @param[in] W Existing wells object. + explicit WellsManager(struct Wells* W); - /// Construct from input deck and grid. - /// The permeability argument may be zero if the input contain - /// well productivity indices, otherwise it must be given in - /// order to approximate these by the Peaceman formula. - WellsManager(const Opm::EclipseGridParser& deck, - const UnstructuredGrid& grid, - const double* permeability); + /// Construct from input deck and grid. + /// The permeability argument may be zero if the input contain + /// well productivity indices, otherwise it must be given in + /// order to approximate these by the Peaceman formula. + WellsManager(const Opm::EclipseStateConstPtr eclipseState, + const size_t timeStep, + const UnstructuredGrid& grid, + const double* permeability); - /// Destructor. - ~WellsManager(); + /// Destructor. + ~WellsManager(); /// Does the "deck" define any wells? bool empty() const; - /// Access the managed Wells. - /// The method is named similarly to c_str() in std::string, - /// to make it clear that we are returning a C-compatible struct. - const Wells* c_wells() const; + /// Access the managed Wells. + /// The method is named similarly to c_str() in std::string, + /// to make it clear that we are returning a C-compatible struct. + const Wells* c_wells() const; /// Access the well group hierarchy. const WellCollection& wellCollection() const; @@ -108,18 +127,31 @@ namespace Opm void applyExplicitReinjectionControls(const std::vector& well_reservoirrates_phase, const std::vector& well_surfacerates_phase); + private: - // Disable copying and assignment. - WellsManager(const WellsManager& other); - WellsManager& operator=(const WellsManager& other); + // Disable copying and assignment. + WellsManager(const WellsManager& other); + WellsManager& operator=(const WellsManager& other); + static void setupCompressedToCartesian(const UnstructuredGrid& grid, std::map& cartesian_to_compressed ); + void setupWellControls(std::vector& wells, size_t timeStep, + std::vector& well_names, const PhaseUsage& phaseUsage); - // Data - Wells* w_; + void 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, + const std::map cartesian_to_compressed, + const double* permeability); + + void addChildGroups(GroupTreeNodeConstPtr parentNode, ScheduleConstPtr schedule, size_t timeStep, const PhaseUsage& phaseUsage); + void setupGuideRates(std::vector& wells, const size_t timeStep, std::vector& well_data, std::map& well_names_to_index); + + + // Data + Wells* w_; WellCollection well_collection_; - - - - }; } // namespace Opm diff --git a/opm/core/wells/well_controls.c b/opm/core/wells/well_controls.c index a1a0b9b00..e4b33a356 100644 --- a/opm/core/wells/well_controls.c +++ b/opm/core/wells/well_controls.c @@ -98,6 +98,8 @@ struct WellControls */ int current; + bool well_is_open; + /* The capacity allocated. */ @@ -130,7 +132,7 @@ well_controls_create(void) ctrl = malloc(1 * sizeof *ctrl); if (ctrl != NULL) { - /* Initialise empty control set */ + /* Initialise empty control set; the well is created open. */ ctrl->num = 0; ctrl->number_of_phases = 0; ctrl->type = NULL; @@ -138,6 +140,7 @@ well_controls_create(void) ctrl->distr = NULL; ctrl->current = -1; ctrl->cpty = 0; + ctrl->well_is_open = true; } return ctrl; @@ -192,11 +195,23 @@ well_controls_set_current( struct WellControls * ctrl, int current) { ctrl->current = current; } -void -well_controls_invert_current( struct WellControls * ctrl ) { - ctrl->current = ~ctrl->current; +bool well_controls_well_is_shut(const struct WellControls * ctrl) { + return !ctrl->well_is_open; } +bool well_controls_well_is_open(const struct WellControls * ctrl) { + return ctrl->well_is_open; +} + +void well_controls_open_well( struct WellControls * ctrl) { + ctrl->well_is_open = true; +} + +void well_controls_shut_well( struct WellControls * ctrl) { + ctrl->well_is_open = false; +} + + enum WellControlType well_controls_iget_type(const struct WellControls * ctrl, int control_index) { @@ -293,19 +308,38 @@ well_controls_add_new(enum WellControlType type , double target , const double * bool -well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2) +well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2 , bool verbose) /* ---------------------------------------------------------------------- */ { bool are_equal = true; - are_equal = (ctrls1->num == ctrls2->num); - are_equal = are_equal && (ctrls1->number_of_phases == ctrls2->number_of_phases); + + if (ctrls1->num != ctrls2->num) { + are_equal = false; + if (verbose) + printf("ctrls1->num:%d ctrls2->num:%d \n",ctrls1->num , ctrls2->num); + } + + if (ctrls1->number_of_phases != ctrls2->number_of_phases) { + are_equal = false; + if (verbose) + printf("ctrls1->number_of_phases:%d ctrls2->number_of_phases:%d \n",ctrls1->number_of_phases , ctrls2->number_of_phases); + } + if (!are_equal) { return are_equal; } - are_equal = are_equal && (memcmp(ctrls1->type, ctrls2->type, ctrls1->num * sizeof *ctrls1->type ) == 0); - are_equal = are_equal && (memcmp(ctrls1->target, ctrls2->target, ctrls1->num * sizeof *ctrls1->target ) == 0); - are_equal = are_equal && (memcmp(ctrls1->distr, ctrls2->distr, ctrls1->num * ctrls1->number_of_phases * sizeof *ctrls1->distr ) == 0); + if (memcmp(ctrls1->type, ctrls2->type, ctrls1->num * sizeof *ctrls1->type ) != 0) { + are_equal = false; + if (verbose) + printf("The ->type vectors are different \n"); + } + + if (memcmp(ctrls1->target, ctrls2->target, ctrls1->num * sizeof *ctrls1->target ) != 0) { + are_equal = false; + if (verbose) + printf("The ->target vectors are different \n"); + } are_equal = are_equal && (ctrls1->cpty == ctrls2->cpty); return are_equal; diff --git a/opm/core/wells/wells.c b/opm/core/wells/wells.c index 13b8c69b3..81a8ddc16 100644 --- a/opm/core/wells/wells.c +++ b/opm/core/wells/wells.c @@ -541,7 +541,7 @@ clone_wells(const struct Wells *W) /* ---------------------------------------------------------------------- */ bool -wells_equal(const struct Wells *W1, const struct Wells *W2) +wells_equal(const struct Wells *W1, const struct Wells *W2 , bool verbose) /* ---------------------------------------------------------------------- */ { bool are_equal = true; @@ -567,10 +567,25 @@ wells_equal(const struct Wells *W1, const struct Wells *W2) are_equal = are_equal && (strcmp(W1->name[i], W2->name[i]) == 0); else are_equal = are_equal && (W1->name[i] == W2->name[i]); + + if (verbose && !are_equal) + printf("Well name[%d] %s and %s are different \n", i , W1->name[i] , W2->name[i]); + } + if (W1->type[i] != W2->type[i]) { + are_equal = false; + if (verbose) + printf("Well->type[%d] different %d %d \n",i , W1->type[i] , W2->type[i] ); + } + if (W1->depth_ref[i] != W2->depth_ref[i]) { + are_equal = false; + if (verbose) + printf("Well->depth_ref[%d] different %g %g \n",i , W1->depth_ref[i] , W2->depth_ref[i] ); + } + if (!well_controls_equal(W1->ctrls[i], W2->ctrls[i],verbose)) { + are_equal = false; + if (verbose) + printf("Well controls are different for well[%d]:%s \n",i,W1->name[i]); } - are_equal = are_equal && (W1->type[i] == W2->type[i]); - are_equal = are_equal && (W1->depth_ref[i] == W2->depth_ref[i]); - are_equal = are_equal && (well_controls_equal(W1->ctrls[i], W2->ctrls[i])); } @@ -581,10 +596,16 @@ wells_equal(const struct Wells *W1, const struct Wells *W2) are_equal = are_equal && (mgmt1->well_cpty == mgmt2->well_cpty); } - are_equal = are_equal && (memcmp(W1->comp_frac, W2->comp_frac, W1->number_of_wells * W1->number_of_phases * sizeof *W1->comp_frac ) == 0); - are_equal = are_equal && (memcmp(W1->well_connpos, W2->well_connpos, (1 + W1->number_of_wells) * sizeof *W1->well_connpos ) == 0); - if (!are_equal) { - return are_equal; + if (memcmp(W1->comp_frac, W2->comp_frac, W1->number_of_wells * W1->number_of_phases * sizeof *W1->comp_frac ) != 0) { + are_equal = false; + if (verbose) + printf("Component fractions different \n"); + } + + if (memcmp(W1->well_connpos, W2->well_connpos, (1 + W1->number_of_wells) * sizeof *W1->well_connpos ) != 0) { + are_equal = false; + if (verbose) + printf("perforation position map difference \n"); } { diff --git a/tests/testBlackoilState1.DATA b/tests/testBlackoilState1.DATA index 5bbb4eb0b..4220b5b45 100644 --- a/tests/testBlackoilState1.DATA +++ b/tests/testBlackoilState1.DATA @@ -10,5 +10,9 @@ DYV DZV 10*10 / +ACTNUM + 1 998*2 3 / + + DEPTHZ 121*2000 / diff --git a/tests/test_linearsolver.cpp b/tests/test_linearsolver.cpp new file mode 100644 index 000000000..184b7765d --- /dev/null +++ b/tests/test_linearsolver.cpp @@ -0,0 +1,186 @@ +/* + Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif +#define NVERBOSE // to suppress our messages when throwing + +#define BOOST_TEST_MODULE OPM-IterativeSolverTest +#include + +#include +#include + +#include +#include +#include +#include + +struct MyMatrix +{ + MyMatrix(int rows, int nnz) + : data(nnz, 0.0), rowStart(rows+1, -1), + colIndex(nnz, -1) + {} + MyMatrix() + : data(), rowStart(), colIndex() + {} + + std::vector data; + std::vector rowStart; + std::vector colIndex; +}; + +std::shared_ptr createLaplacian(int N) +{ + MyMatrix* mm=new MyMatrix(N*N, N*N*5); + int nnz=0; + mm->rowStart[0]=0; + for(int row=0; row0) + { + mm->colIndex[nnz]=row-N; + mm->data[nnz++]=-1; + } + if(x>0) + { + mm->colIndex[nnz]=row-1; + mm->data[nnz++]=-1; + } + mm->colIndex[nnz]=row; + mm->data[nnz++]=4; + if(xcolIndex[nnz]=row+1; + mm->data[nnz++]=-1; + } + if(ycolIndex[nnz]=row+N; + mm->data[nnz++]=-1; + } + mm->rowStart[row+1]=nnz; + } + mm->data.resize(nnz); + mm->colIndex.resize(nnz); + return std::shared_ptr(mm); +} + +void createRandomVectors(int NN, std::vector& x, std::vector& b, + const MyMatrix& mat) +{ + x.resize(NN); + for(auto entry=x.begin(), end =x.end(); entry!=end; ++entry) + *entry=((double) (rand()%100))/10.0; + + b.resize(NN); + std::fill(b.begin(), b.end(), 0.0); + + // Construct the right hand side as b=A*x + for(std::size_t row=0; row x(N*N), b(N*N); + createRandomVectors(100*100, x, b, *mat); + std::vector exact(x); + std::fill(x.begin(), x.end(), 0.0); + Opm::LinearSolverFactory ls(param); + ls.solve(N*N, mat->data.size(), &(mat->rowStart[0]), + &(mat->colIndex[0]), &(mat->data[0]), &(b[0]), + &(x[0])); +} + + +BOOST_AUTO_TEST_CASE(DefaultTest) +{ + Opm::parameter::ParameterGroup param; + param.insertParameter(std::string("linsolver_max_iterations"), std::string("200")); + param.insertParameter(std::string("linsolver_verbosity"), std::string("2")); + run_test(param); +} + +#ifdef HAVE_DUNE_ISTL +BOOST_AUTO_TEST_CASE(CGAMGTest) +{ + Opm::parameter::ParameterGroup param; + param.insertParameter(std::string("linsolver"), std::string("istl")); + param.insertParameter(std::string("linsolver_type"), std::string("1")); + param.insertParameter(std::string("linsolver_max_iterations"), std::string("200")); + param.insertParameter(std::string("linsolver_verbosity"), std::string("2")); + run_test(param); +} + +BOOST_AUTO_TEST_CASE(CGILUTest) +{ + Opm::parameter::ParameterGroup param; + param.insertParameter(std::string("linsolver"), std::string("istl")); + param.insertParameter(std::string("linsolver_type"), std::string("0")); + param.insertParameter(std::string("linsolver_max_iterations"), std::string("200")); + param.insertParameter(std::string("linsolver_verbosity"), std::string("2")); + run_test(param); +} + +BOOST_AUTO_TEST_CASE(BiCGILUTest) +{ + Opm::parameter::ParameterGroup param; + param.insertParameter(std::string("linsolver"), std::string("istl")); + param.insertParameter(std::string("linsolver_type"), std::string("2")); + param.insertParameter(std::string("linsolver_max_iterations"), std::string("200")); + param.insertParameter(std::string("linsolver_verbosity"), std::string("2")); + run_test(param); +} + +#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) +BOOST_AUTO_TEST_CASE(FastAMGTest) +{ + Opm::parameter::ParameterGroup param; + param.insertParameter(std::string("linsolver"), std::string("istl")); + param.insertParameter(std::string("linsolver_type"), std::string("3")); + param.insertParameter(std::string("linsolver_max_iterations"), std::string("200")); + param.insertParameter(std::string("linsolver_verbosity"), std::string("2")); + run_test(param); +} + +BOOST_AUTO_TEST_CASE(KAMGTest) +{ + Opm::parameter::ParameterGroup param; + param.insertParameter(std::string("linsolver"), std::string("istl")); + param.insertParameter(std::string("linsolver_type"), std::string("4")); + param.insertParameter(std::string("linsolver_max_iterations"), std::string("200")); + run_test(param); +} +#endif +#endif diff --git a/tests/test_wellcollection.cpp b/tests/test_wellcollection.cpp new file mode 100644 index 000000000..d22ba3b4b --- /dev/null +++ b/tests/test_wellcollection.cpp @@ -0,0 +1,91 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE WellCollectionTest +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace Opm; + +BOOST_AUTO_TEST_CASE(AddWellsAndGroupToCollection) { + ParserPtr parser(new Parser()); + boost::filesystem::path scheduleFile("wells_group.data"); + DeckConstPtr deck = parser->parseFile(scheduleFile.string()); + EclipseStateConstPtr eclipseState(new EclipseState(deck)); + PhaseUsage pu = phaseUsageFromDeck(eclipseState); + + GroupTreeNodePtr field=eclipseState->getSchedule()->getGroupTree(2)->getNode("FIELD"); + GroupTreeNodePtr g1=eclipseState->getSchedule()->getGroupTree(2)->getNode("G1"); + GroupTreeNodePtr g2=eclipseState->getSchedule()->getGroupTree(2)->getNode("G2"); + + WellCollection collection; + + // Add groups to WellCollection + GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(field->name()); + collection.addField(fieldGroup, 2, pu); + + for (auto iter = field->begin(); iter != field->end(); ++iter) { + GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name()); + collection.addGroup(childGroupNode, fieldGroup->name(), 2, pu); + } + + GroupConstPtr g1Group = eclipseState->getSchedule()->getGroup(g1->name()); + for (auto iter = g1->begin(); iter != g1->end(); ++iter) { + GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name()); + collection.addGroup(childGroupNode, g1Group->name(), 2, pu); + } + + + GroupConstPtr g2Group = eclipseState->getSchedule()->getGroup(g2->name()); + for (auto iter = g2->begin(); iter != g2->end(); ++iter) { + GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name()); + collection.addGroup(childGroupNode, g2Group->name(), 2, pu); + } + + BOOST_CHECK_EQUAL("FIELD", collection.findNode("FIELD")->name()); + BOOST_CHECK_EQUAL("FIELD", collection.findNode("G1")->getParent()->name()); + BOOST_CHECK_EQUAL("FIELD", collection.findNode("G2")->getParent()->name()); + + // Add wells to WellCollection + WellCollection wellCollection; + std::vector wells = eclipseState->getSchedule()->getWells(); + for (size_t i=0; igetParent()->name()); + BOOST_CHECK_EQUAL("G1", collection.findNode("INJ2")->getParent()->name()); + BOOST_CHECK_EQUAL("G2", collection.findNode("PROD1")->getParent()->name()); + BOOST_CHECK_EQUAL("G2", collection.findNode("PROD2")->getParent()->name()); +} + diff --git a/tests/test_wellcontrols.cpp b/tests/test_wellcontrols.cpp index 096b9c6a0..8424ab0b5 100644 --- a/tests/test_wellcontrols.cpp +++ b/tests/test_wellcontrols.cpp @@ -45,11 +45,6 @@ BOOST_AUTO_TEST_CASE(Construction) well_controls_set_current( ctrls , 2 ); BOOST_CHECK_EQUAL( 2 , well_controls_get_current( ctrls )); - well_controls_invert_current( ctrls ); - BOOST_CHECK( well_controls_get_current( ctrls ) < 0 ); - well_controls_invert_current( ctrls ); - BOOST_CHECK_EQUAL( 2 , well_controls_get_current( ctrls )); - { enum WellControlType type1 = BHP; enum WellControlType type2 = SURFACE_RATE; @@ -103,3 +98,21 @@ BOOST_AUTO_TEST_CASE(Construction) } +BOOST_AUTO_TEST_CASE(OpenClose) +{ + struct WellControls * ctrls = well_controls_create(); + + BOOST_CHECK_EQUAL( true , well_controls_well_is_open(ctrls) ); + BOOST_CHECK_EQUAL( false , well_controls_well_is_shut(ctrls) ); + + well_controls_open_well( ctrls ); + BOOST_CHECK_EQUAL( true , well_controls_well_is_open(ctrls) ); + BOOST_CHECK_EQUAL( false , well_controls_well_is_shut(ctrls) ); + + well_controls_shut_well( ctrls ); + BOOST_CHECK_EQUAL( false , well_controls_well_is_open(ctrls) ); + BOOST_CHECK_EQUAL( true , well_controls_well_is_shut(ctrls) ); + + well_controls_destroy( ctrls ); +} + diff --git a/tests/test_wells.cpp b/tests/test_wells.cpp index 9e76d899d..8502c4b87 100644 --- a/tests/test_wells.cpp +++ b/tests/test_wells.cpp @@ -204,7 +204,7 @@ BOOST_AUTO_TEST_CASE(Copy) BOOST_CHECK_EQUAL( dist1[p] , dist2[p]); } } - BOOST_CHECK( well_controls_equal( c1 , c2 )); + BOOST_CHECK( well_controls_equal( c1 , c2 , false) ); } } } @@ -219,7 +219,7 @@ BOOST_AUTO_TEST_CASE(Equals_WellsEqual_ReturnsTrue) { std::shared_ptr W2(create_wells(nphases, nwells, nperfs), destroy_wells); - BOOST_CHECK(wells_equal(W1.get(), W2.get())); + BOOST_CHECK(wells_equal(W1.get(), W2.get() , false)); } @@ -232,5 +232,5 @@ BOOST_AUTO_TEST_CASE(Equals_WellsDiffer_ReturnsFalse) { std::shared_ptr W2(create_wells(nphases, 3, nperfs), destroy_wells); - BOOST_CHECK(!wells_equal(W1.get(), W2.get())); + BOOST_CHECK(!wells_equal(W1.get(), W2.get() , false )); } diff --git a/tests/test_wellsgroup.cpp b/tests/test_wellsgroup.cpp new file mode 100644 index 000000000..46fd8eb0c --- /dev/null +++ b/tests/test_wellsgroup.cpp @@ -0,0 +1,110 @@ +/* + Copyright 2014 Statoil. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE WellsGroupTest + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +using namespace Opm; + +BOOST_AUTO_TEST_CASE(ConstructGroupFromWell) { + ParserPtr parser(new Parser()); + boost::filesystem::path scheduleFile("wells_group.data"); + DeckConstPtr deck = parser->parseFile(scheduleFile.string()); + EclipseStateConstPtr eclipseState(new EclipseState(deck)); + PhaseUsage pu = phaseUsageFromDeck(eclipseState); + + std::vector wells = eclipseState->getSchedule()->getWells(); + + for (size_t i=0; i wellsGroup = createWellWellsGroup(well, 2, pu); + BOOST_CHECK_EQUAL(well->name(), wellsGroup->name()); + if (well->isInjector(2)) { + const WellInjectionProperties& properties = well->getInjectionProperties(2); + BOOST_CHECK_EQUAL(properties.surfaceInjectionRate, wellsGroup->injSpec().surface_flow_max_rate_); + BOOST_CHECK_EQUAL(properties.BHPLimit, wellsGroup->injSpec().BHP_limit_); + BOOST_CHECK_EQUAL(properties.reservoirInjectionRate, wellsGroup->injSpec().reservoir_flow_max_rate_); + BOOST_CHECK_EQUAL(0.0, wellsGroup->prodSpec().guide_rate_); + } + if (well->isProducer(2)) { + const WellProductionProperties& properties = well->getProductionProperties(2); + BOOST_CHECK_EQUAL(properties.ResVRate, wellsGroup->prodSpec().reservoir_flow_max_rate_); + BOOST_CHECK_EQUAL(properties.BHPLimit, wellsGroup->prodSpec().BHP_limit_); + BOOST_CHECK_EQUAL(properties.OilRate, wellsGroup->prodSpec().oil_max_rate_); + BOOST_CHECK_EQUAL(properties.WaterRate, wellsGroup->prodSpec().water_max_rate_); + BOOST_CHECK_EQUAL(0.0, wellsGroup->injSpec().guide_rate_); + } + } +} + + +BOOST_AUTO_TEST_CASE(ConstructGroupFromGroup) { + ParserPtr parser(new Parser()); + boost::filesystem::path scheduleFile("wells_group.data"); + DeckConstPtr deck = parser->parseFile(scheduleFile.string()); + EclipseStateConstPtr eclipseState(new EclipseState(deck)); + PhaseUsage pu = phaseUsageFromDeck(eclipseState); + + std::vector nodes = eclipseState->getSchedule()->getGroupTree(2)->getNodes(); + + for (size_t i=0; igetSchedule()->getGroup(nodes[i]->name()); + std::shared_ptr wellsGroup = createGroupWellsGroup(group, 2, pu); + BOOST_CHECK_EQUAL(group->name(), wellsGroup->name()); + if (group->isInjectionGroup(2)) { + BOOST_CHECK_EQUAL(group->getSurfaceMaxRate(2), wellsGroup->injSpec().surface_flow_max_rate_); + BOOST_CHECK_EQUAL(group->getReservoirMaxRate(2), wellsGroup->injSpec().reservoir_flow_max_rate_); + BOOST_CHECK_EQUAL(group->getTargetReinjectFraction(2), wellsGroup->injSpec().reinjection_fraction_target_); + BOOST_CHECK_EQUAL(group->getTargetVoidReplacementFraction(2), wellsGroup->injSpec().voidage_replacment_fraction_); + } + if (group->isProductionGroup(2)) { + BOOST_CHECK_EQUAL(group->getReservoirMaxRate(2), wellsGroup->prodSpec().reservoir_flow_max_rate_); + BOOST_CHECK_EQUAL(group->getGasTargetRate(2), wellsGroup->prodSpec().gas_max_rate_); + BOOST_CHECK_EQUAL(group->getOilTargetRate(2), wellsGroup->prodSpec().oil_max_rate_); + BOOST_CHECK_EQUAL(group->getWaterTargetRate(2), wellsGroup->prodSpec().water_max_rate_); + BOOST_CHECK_EQUAL(group->getLiquidTargetRate(2), wellsGroup->prodSpec().liquid_max_rate_); + } + } +} + + + diff --git a/tests/test_wellsmanager.cpp b/tests/test_wellsmanager.cpp index b3465fac6..2141e4bf8 100644 --- a/tests/test_wellsmanager.cpp +++ b/tests/test_wellsmanager.cpp @@ -27,6 +27,10 @@ #define BOOST_TEST_MODULE WellsManagerTests #include + +#include +#include + #include #include #include @@ -57,7 +61,7 @@ void wells_static_check(const Wells* wells) { } -/* +/* The number of controls is determined by looking at which elements have been given explicit - non-default - values in the WCONxxxx keyword. Is that at all interesting? @@ -67,7 +71,7 @@ void wells_static_check(const Wells* wells) { void check_controls_epoch0( struct WellControls ** ctrls) { // The injector { - const struct WellControls * ctrls0 = ctrls[0]; + const struct WellControls * ctrls0 = ctrls[0]; BOOST_CHECK_EQUAL( 3 , well_controls_get_num(ctrls0)); // The number of controls for the injector == 3?? BOOST_CHECK_EQUAL( SURFACE_RATE , well_controls_iget_type(ctrls0 , 0) ); @@ -83,7 +87,7 @@ void check_controls_epoch0( struct WellControls ** ctrls) { BOOST_CHECK_EQUAL( 0 , well_controls_get_current(ctrls0) ); // The phase distribution in the active target - { + { const double * distr = well_controls_iget_distr( ctrls0 , 0 ); BOOST_CHECK_EQUAL( 0 , distr[0] ); // Water BOOST_CHECK_EQUAL( 0 , distr[1] ); // Oil @@ -106,7 +110,7 @@ void check_controls_epoch0( struct WellControls ** ctrls) { BOOST_CHECK_EQUAL( 0 , well_controls_get_current(ctrls1)); // The phase distribution in the active target - { + { const double * distr = well_controls_iget_distr( ctrls1 , 0 ); BOOST_CHECK_EQUAL( 0 , distr[0] ); // Water BOOST_CHECK_EQUAL( 1 , distr[1] ); // Oil @@ -121,7 +125,7 @@ void check_controls_epoch0( struct WellControls ** ctrls) { void check_controls_epoch1( struct WellControls ** ctrls) { // The injector { - const struct WellControls * ctrls0 = ctrls[0]; + const struct WellControls * ctrls0 = ctrls[0]; BOOST_CHECK_EQUAL( 3 , well_controls_get_num(ctrls0)); // The number of controls for the injector == 3?? BOOST_CHECK_EQUAL( SURFACE_RATE , well_controls_iget_type(ctrls0 , 0 )); @@ -136,7 +140,7 @@ void check_controls_epoch1( struct WellControls ** ctrls) { // Which control is active BOOST_CHECK_EQUAL( 1 , well_controls_get_current(ctrls0)); - { + { const double * distr = well_controls_iget_distr( ctrls0 , 1 ); BOOST_CHECK_EQUAL( 1 , distr[0] ); // Water BOOST_CHECK_EQUAL( 0 , distr[1] ); // Oil @@ -160,7 +164,7 @@ void check_controls_epoch1( struct WellControls ** ctrls) { // Which control is active BOOST_CHECK_EQUAL( 1 , well_controls_get_current(ctrls1) ); - { + { const double * distr = well_controls_iget_distr( ctrls1 , 1 ); BOOST_CHECK_EQUAL( 1 , distr[0] ); // Water BOOST_CHECK_EQUAL( 1 , distr[1] ); // Oil @@ -169,29 +173,24 @@ void check_controls_epoch1( struct WellControls ** ctrls) { } } +BOOST_AUTO_TEST_CASE(New_Constructor_Works) { + Opm::ParserPtr parser(new Opm::Parser()); + Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data.data"))); - -BOOST_AUTO_TEST_CASE(Constructor_Works) { Opm::EclipseGridParser Deck("wells_manager_data.data"); Opm::GridManager gridManager(Deck); - Deck.setCurrentEpoch(0); { - Opm::WellsManager wellsManager(Deck, *gridManager.c_grid(), NULL); - const Wells* wells = wellsManager.c_wells(); - wells_static_check( wells ); - check_controls_epoch0( wells->ctrls ); + Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL); + wells_static_check( wellsManager.c_wells() ); + check_controls_epoch0( wellsManager.c_wells()->ctrls ); } - - Deck.setCurrentEpoch(1); { - Opm::WellsManager wellsManager(Deck, *gridManager.c_grid(), NULL); - const Wells* wells = wellsManager.c_wells(); - - wells_static_check( wells ); - check_controls_epoch1( wells->ctrls ); + Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL); + wells_static_check( wellsManager.c_wells() ); + check_controls_epoch1( wellsManager.c_wells()->ctrls ); } } @@ -200,15 +199,14 @@ BOOST_AUTO_TEST_CASE(Constructor_Works) { BOOST_AUTO_TEST_CASE(WellsEqual) { Opm::EclipseGridParser Deck("wells_manager_data.data"); Opm::GridManager gridManager(Deck); + Opm::ParserPtr parser(new Opm::Parser()); + Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data.data"))); - Deck.setCurrentEpoch(0); - Opm::WellsManager wellsManager0(Deck, *gridManager.c_grid(), NULL); + Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL); + Opm::WellsManager wellsManager1(eclipseState , 1 , *gridManager.c_grid(), NULL); - Deck.setCurrentEpoch(1); - Opm::WellsManager wellsManager1(Deck, *gridManager.c_grid(), NULL); - - BOOST_CHECK( wells_equal( wellsManager0.c_wells() , wellsManager0.c_wells()) ); - BOOST_CHECK( !wells_equal( wellsManager0.c_wells() , wellsManager1.c_wells()) ); + BOOST_CHECK( wells_equal( wellsManager0.c_wells() , wellsManager0.c_wells(),false) ); + BOOST_CHECK( !wells_equal( wellsManager0.c_wells() , wellsManager1.c_wells(),false) ); } @@ -216,21 +214,36 @@ BOOST_AUTO_TEST_CASE(ControlsEqual) { Opm::EclipseGridParser Deck("wells_manager_data.data"); Opm::GridManager gridManager(Deck); - Deck.setCurrentEpoch(0); - Opm::WellsManager wellsManager0(Deck, *gridManager.c_grid(), NULL); + Opm::ParserPtr parser(new Opm::Parser()); + Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data.data"))); - Deck.setCurrentEpoch(1); - Opm::WellsManager wellsManager1(Deck, *gridManager.c_grid(), NULL); + Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL); + Opm::WellsManager wellsManager1(eclipseState , 1 , *gridManager.c_grid(), NULL); - BOOST_CHECK( well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0])); - BOOST_CHECK( well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1])); - BOOST_CHECK( well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager1.c_wells()->ctrls[0])); - BOOST_CHECK( well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager1.c_wells()->ctrls[1])); + BOOST_CHECK( well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false)); + BOOST_CHECK( well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false)); + BOOST_CHECK( well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager1.c_wells()->ctrls[0] , false)); + BOOST_CHECK( well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager1.c_wells()->ctrls[1] , false)); - BOOST_CHECK( !well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[1])); - BOOST_CHECK( !well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[0])); - BOOST_CHECK( !well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0])); - BOOST_CHECK( !well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1])); + BOOST_CHECK( !well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[1] , false)); + BOOST_CHECK( !well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[0] , false)); + BOOST_CHECK( !well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false)); + BOOST_CHECK( !well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false)); } + +BOOST_AUTO_TEST_CASE(WellHasSTOP_ExceptionIsThrown) { + Opm::EclipseGridParser Deck("wells_manager_data_wellSTOP.data"); + Opm::GridManager gridManager(Deck); + + Opm::ParserPtr parser(new Opm::Parser()); + Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data_wellSTOP.data"))); + + Deck.setCurrentEpoch(0); + + BOOST_CHECK_THROW( new Opm::WellsManager(eclipseState, 0, *gridManager.c_grid(), NULL), std::runtime_error ); +} + + + diff --git a/tests/wells_group.data b/tests/wells_group.data new file mode 100755 index 000000000..5d2e58cb2 --- /dev/null +++ b/tests/wells_group.data @@ -0,0 +1,68 @@ +OIL +GAS +WATER + +DIMENS + 10 10 5 / + +GRID + +DXV +10*1000.0 / + +DYV +10*1000.0 / + +DZV +10.0 20.0 30.0 10.0 5.0 / + +DEPTHZ +121*2000 +/ + +SCHEDULE + +GRUPTREE + 'G1' 'FIELD' / + 'G2' 'FIELD' / +/ + + +WELSPECS + 'INJ1' 'G1' 1 1 8335 'GAS' / + 'PROD1' 'G2' 10 10 8400 'OIL' / +/ + +TSTEP + 14.0 / +/ + +WELSPECS + 'INJ2' 'G1' 1 1 8335 'GAS' / + 'PROD2' 'G2' 10 10 8400 'OIL' / +/ + +GCONINJE +'G1' GAS RATE 30000 / +/ + +GCONPROD +'G2' ORAT 10000 / +/ + +WCONINJE + 'INJ1' 'WATER' 'OPEN' 'RESV' 10 20 40 / + 'INJ2' 'WATER' 'OPEN' 'RESV' 10 20 40 / +/ + +WCONPROD + 'PROD1' 'OPEN' 'RESV' 999 3* 123 100 / + 'PROD2' 'OPEN' 'RESV' 999 3* 123 100 / +/ + +TSTEP + 3 / +/ + + +END diff --git a/tests/wells_manager_data.data b/tests/wells_manager_data.data index f22152c13..732b528be 100755 --- a/tests/wells_manager_data.data +++ b/tests/wells_manager_data.data @@ -17,8 +17,8 @@ DZV 10.0 20.0 30.0 10.0 5.0 / DEPTHZ -121*2000 -/ +121*2000 / + SCHEDULE @@ -34,11 +34,11 @@ COMPDAT WCONPROD 'PROD1' 'OPEN' 'ORAT' 20000 4* 1000 / - / +/ WCONINJE 'INJ1' 'GAS' 'OPEN' 'RATE' 100 200 400 / - / +/ DATES @@ -53,10 +53,26 @@ WCONINJE 'INJ1' 'WATER' 'OPEN' 'RESV' 10 20 40 / / - +END TSTEP -14.0 - / + 14.0 / +/ + + +WELSPECS + 'TEST1' 'G1' 1 1 8335 'GAS' / + 'TEST2' 'G2' 10 10 8400 'OIL' / +/ + +GRUPTREE + 'G1' 'SHIP' / + 'G2' 'SHIP' / +/ + +TSTEP + 3 / +/ + END diff --git a/tests/wells_manager_data_expanded.data b/tests/wells_manager_data_expanded.data new file mode 100755 index 000000000..672077681 --- /dev/null +++ b/tests/wells_manager_data_expanded.data @@ -0,0 +1,86 @@ +OIL +GAS +WATER + +DIMENS + 10 10 5 / + +GRID + +DXV +10*1000.0 / + +DYV +10*1000.0 / + +DZV +10.0 20.0 30.0 10.0 5.0 / + +-- The DEPTHZ keyword is only here to satisfy the old parser; content might +-- completely bogus. +DEPTHZ +121*2000 / + + +SCHEDULE + +WELSPECS + 'INJ1' 'G' 1 1 8335 'GAS' / + 'PROD1' 'G' 10 10 8400 'OIL' / +/ + +COMPDAT + 'INJ1' 1 1 1 2 'OPEN' 1 10.6092 0.5 / + 'INJ1' 1 1 3 5 'OPEN' 1 12.6092 0.5 / + 'INJ1' 2 2 1 1 'OPEN' 1 14.6092 0.5 / + 'PROD1' 10 3 3 3 'OPEN' 0 10.6092 0.5 / +/ + +WCONPROD + 'PROD1' 'OPEN' 'ORAT' 20000 4* 1000 / +/ + +WCONINJE + 'INJ1' 'GAS' 'OPEN' 'RATE' 100 200 400 / +/ + + +DATES + 1 'FEB' 2000 / +/ + +WELSPECS + 'INJ1' 'G3' 1 1 8335 'GAS' / + 'QNJ2' 'G3' 1 1 8335 'GAS' / +/ + + +COMPDAT + 'QNJ2' 3 4 1 2 'OPEN' 1 10.6092 0.5 / + 'QNJ2' 4 4 3 5 'OPEN' 1 12.6092 0.5 / +/ + +WCONPROD + 'PROD1' 'OPEN' 'RESV' 999 3* 123 100 / +/ + +WCONINJE + 'INJ1' 'WATER' 'OPEN' 'RESV' 10 20 40 / + 'QNJ2' 'WATER' 'OPEN' 'RESV' 7 33 39 / +/ + + + +TSTEP +14.0 / +/ + +WELOPEN + 'INJ1' 'SHUT' 5* / +/ + +TSTEP +14.0 / +/ + +END diff --git a/tests/wells_manager_data_wellSTOP.data b/tests/wells_manager_data_wellSTOP.data new file mode 100755 index 000000000..bdcbc3817 --- /dev/null +++ b/tests/wells_manager_data_wellSTOP.data @@ -0,0 +1,42 @@ +OIL +GAS +WATER + +DIMENS + 10 10 5 / + +GRID + +DXV +10*1000.0 / + +DYV +10*1000.0 / + +DZV +10.0 20.0 30.0 10.0 5.0 / + +DEPTHZ +121*2000 +/ + +SCHEDULE + +WELSPECS + 'INJ1' 'G' 1 1 8335 'GAS' / + 'PROD1' 'G' 10 10 8400 'OIL' / +/ + +COMPDAT + 'INJ1' 1 1 1 1 'OPEN' 1 10.6092 0.5 / + 'PROD1' 10 3 3 3 'OPEN' 0 10.6092 0.5 / +/ + +WELOPEN + 'INJ1' 'STOP' 5* / +/ + +TSTEP + 10 / + +END