From 7fdda982136e125e271cab571838ec28938907d7 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 24 Jan 2014 12:04:04 +0100 Subject: [PATCH 1/4] add variants of all methods which take a deck of the new parser to the rock properties --- opm/core/props/rock/RockCompressibility.cpp | 49 ++++ opm/core/props/rock/RockCompressibility.hpp | 6 + opm/core/props/rock/RockFromDeck.cpp | 250 ++++++++++++++++++++ opm/core/props/rock/RockFromDeck.hpp | 16 ++ 4 files changed, 321 insertions(+) diff --git a/opm/core/props/rock/RockCompressibility.cpp b/opm/core/props/rock/RockCompressibility.cpp index 7b0f69cb8..5b990af14 100644 --- a/opm/core/props/rock/RockCompressibility.cpp +++ b/opm/core/props/rock/RockCompressibility.cpp @@ -25,6 +25,8 @@ #include #include +#include + #include namespace Opm @@ -65,6 +67,53 @@ 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."); + + std::vector rtColumnNames + { + "PRESSURE", + "POREVOL MULT", + "TRANSMISC MULT" + }; + if (newParserDeck->hasKeyword("RKTRMDIR")) + { + // the number of colums of the "ROCKTAB" keyword + // depends on the presence of the "RKTRMDIR" + // keyword. Messy stuff... + rtColumnNames.push_back("TRANSMISC MULT Y"); + rtColumnNames.push_back("TRANSMISC MULT Z"); + + // well, okay. we don't support non-isotropic + // transmiscibility multipliers yet + OPM_THROW(std::runtime_error, "Support for non-isotropic " + "transmiscibility multipliers is not implemented yet."); + }; + + Opm::SimpleTable rtTable(rtKeyword, rtColumnNames); + + p_ = rtTable.getColumn(0); + poromult_ = rtTable.getColumn(1); + transmult_ = rtTable.getColumn(2); + } else if (newParserDeck->hasKeyword("ROCK")) { + Opm::DeckKeywordConstPtr rockKeyword = newParserDeck->getKeyword("ROCK"); + if (rockKeyword->size() != 1) + OPM_THROW(std::runtime_error, "Can only handle a single region in ROCK."); + + Opm::DeckRecordConstPtr rockRecord = rockKeyword->getRecord(0); + pref_ = rockRecord->getItem(0)->getSIDouble(0); + rock_comp_ = rockRecord->getItem(1)->getSIDouble(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_; From e177b657f9ed62666f56af02bde9817892fb53c5 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 5 Feb 2014 14:54:13 +0100 Subject: [PATCH 2/4] use the (new) RocktabTable class for rock compressibility thanks to @bska for pointing this out --- opm/core/props/rock/RockCompressibility.cpp | 32 ++++++++------------- 1 file changed, 12 insertions(+), 20 deletions(-) diff --git a/opm/core/props/rock/RockCompressibility.cpp b/opm/core/props/rock/RockCompressibility.cpp index 5b990af14..0f2cefecb 100644 --- a/opm/core/props/rock/RockCompressibility.cpp +++ b/opm/core/props/rock/RockCompressibility.cpp @@ -25,7 +25,7 @@ #include #include -#include +#include #include @@ -76,31 +76,23 @@ namespace Opm if (rtKeyword->size() != 1) OPM_THROW(std::runtime_error, "Can only handle a single region in ROCKTAB."); - std::vector rtColumnNames + // 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) { - "PRESSURE", - "POREVOL MULT", - "TRANSMISC MULT" - }; - if (newParserDeck->hasKeyword("RKTRMDIR")) - { - // the number of colums of the "ROCKTAB" keyword - // depends on the presence of the "RKTRMDIR" - // keyword. Messy stuff... - rtColumnNames.push_back("TRANSMISC MULT Y"); - rtColumnNames.push_back("TRANSMISC MULT Z"); - // well, okay. we don't support non-isotropic - // transmiscibility multipliers yet + // transmisibility multipliers yet OPM_THROW(std::runtime_error, "Support for non-isotropic " - "transmiscibility multipliers is not implemented yet."); + "transmisibility multipliers is not implemented yet."); }; - Opm::SimpleTable rtTable(rtKeyword, rtColumnNames); + Opm::RocktabTable rocktabTable(rtKeyword, isDirectional); - p_ = rtTable.getColumn(0); - poromult_ = rtTable.getColumn(1); - transmult_ = rtTable.getColumn(2); + p_ = rocktabTable.getPressureColumn(); + poromult_ = rocktabTable.getPoreVolumeMultiplierColumn(); + transmult_ = rocktabTable.getTransmisibilityMultiplierColumn(); } else if (newParserDeck->hasKeyword("ROCK")) { Opm::DeckKeywordConstPtr rockKeyword = newParserDeck->getKeyword("ROCK"); if (rockKeyword->size() != 1) From 5b77bcd4a0c10bd547a7ba49a8061cfc61fcc794 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 5 Feb 2014 15:23:05 +0100 Subject: [PATCH 3/4] handle the ROCK keywords using the new opm-parser utility class --- opm/core/props/rock/RockCompressibility.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/opm/core/props/rock/RockCompressibility.cpp b/opm/core/props/rock/RockCompressibility.cpp index 0f2cefecb..1144d701e 100644 --- a/opm/core/props/rock/RockCompressibility.cpp +++ b/opm/core/props/rock/RockCompressibility.cpp @@ -26,6 +26,7 @@ #include #include +#include #include @@ -94,13 +95,12 @@ namespace Opm poromult_ = rocktabTable.getPoreVolumeMultiplierColumn(); transmult_ = rocktabTable.getTransmisibilityMultiplierColumn(); } else if (newParserDeck->hasKeyword("ROCK")) { - Opm::DeckKeywordConstPtr rockKeyword = newParserDeck->getKeyword("ROCK"); - if (rockKeyword->size() != 1) + Opm::RockTable rockTable(newParserDeck->getKeyword("ROCK")); + if (rockTable.numRows() != 1) OPM_THROW(std::runtime_error, "Can only handle a single region in ROCK."); - Opm::DeckRecordConstPtr rockRecord = rockKeyword->getRecord(0); - pref_ = rockRecord->getItem(0)->getSIDouble(0); - rock_comp_ = rockRecord->getItem(1)->getSIDouble(0); + 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; } From 1221ea25c10efe4ee1ec606d7301f4be0c11e557 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 5 Feb 2014 15:58:08 +0100 Subject: [PATCH 4/4] more typos thanks to @bska --- opm/core/props/rock/RockCompressibility.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/core/props/rock/RockCompressibility.cpp b/opm/core/props/rock/RockCompressibility.cpp index 1144d701e..d1729533c 100644 --- a/opm/core/props/rock/RockCompressibility.cpp +++ b/opm/core/props/rock/RockCompressibility.cpp @@ -84,16 +84,16 @@ namespace Opm if (isDirectional) { // well, okay. we don't support non-isotropic - // transmisibility multipliers yet + // transmissibility multipliers yet OPM_THROW(std::runtime_error, "Support for non-isotropic " - "transmisibility multipliers is not implemented yet."); + "transmissibility multipliers is not implemented yet."); }; Opm::RocktabTable rocktabTable(rtKeyword, isDirectional); p_ = rocktabTable.getPressureColumn(); poromult_ = rocktabTable.getPoreVolumeMultiplierColumn(); - transmult_ = rocktabTable.getTransmisibilityMultiplierColumn(); + transmult_ = rocktabTable.getTransmissibilityMultiplierColumn(); } else if (newParserDeck->hasKeyword("ROCK")) { Opm::RockTable rockTable(newParserDeck->getKeyword("ROCK")); if (rockTable.numRows() != 1)