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_;