mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge pull request #467 from andlaus/parser-integrate_rock_properties
add variants of all methods which take a deck of the new parser to the rock properties
This commit is contained in:
Joakim Hove
commit
0393a8f87d
@ -25,6 +25,9 @@
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/linearInterpolation.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Utility/RocktabTable.hpp>
|
||||
#include <opm/parser/eclipse/Utility/RockTable.hpp>
|
||||
|
||||
#include <iostream>
|
||||
|
||||
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);
|
||||
|
||||
@ -20,6 +20,8 @@
|
||||
#ifndef OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED
|
||||
#define OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED
|
||||
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
|
||||
#include <vector>
|
||||
|
||||
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]
|
||||
|
||||
@ -21,6 +21,9 @@
|
||||
#include "config.h"
|
||||
#include <opm/core/props/rock/RockFromDeck.hpp>
|
||||
#include <opm/core/grid.h>
|
||||
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
|
||||
#include <array>
|
||||
|
||||
namespace Opm
|
||||
@ -37,6 +40,10 @@ namespace Opm
|
||||
PermeabilityKind fillTensor(const EclipseGridParser& parser,
|
||||
std::vector<const std::vector<double>*>& tensor,
|
||||
std::array<int,9>& kmap);
|
||||
PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck,
|
||||
std::vector<const std::vector<double>*>& tensor,
|
||||
std::array<int,9>& 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<double>& 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<const std::vector<double>*> tensor;
|
||||
tensor.reserve(10);
|
||||
|
||||
const std::vector<double> zero(num_global_cells, 0.0);
|
||||
tensor.push_back(&zero);
|
||||
|
||||
std::array<int,9> 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<unsigned char>::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<const std::vector<double>*>& tensor,
|
||||
std::array<int,9>& 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
|
||||
|
||||
@ -22,6 +22,9 @@
|
||||
|
||||
|
||||
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
|
||||
#include <vector>
|
||||
|
||||
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<double> porosity_;
|
||||
std::vector<double> permeability_;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user