OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Reimplement assign* in terms of Compressed<>

Browse Source 
This is a demonstration of using the

    GridPropertyAccess::Compressed<>

class template.  We save (some) memory by not creating the zero
fall-back vector in assignPermeability(), preferring instead to use
the fall-back/default mechanism of ArrayPolicy::ExtractFromDeck<>.

While here, adjust vector<PermComponent>::reserve() capacity to
reflect actual requirements.
This commit is contained in:
Bård Skaflestad 2014-08-27 21:14:31 +02:00
parent 37dfc4a3b5
commit 58f55f3c90
1 changed files with 69 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

109
opm/core/props/rock/RockFromDeck.cpp
View File

@ -26,7 +26,11 @@
#include <opm/parser/eclipse/Deck/Deck.hpp> #include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp> #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/core/utility/CompressedPropertyAccess.hpp>
#include <array> #include <array>
#include <string>
#include <vector>
namespace Opm namespace Opm
{ {
@ -38,9 +42,23 @@ namespace Opm
void setScalarPermIfNeeded(std::array<int,9>& kmap, void setScalarPermIfNeeded(std::array<int,9>& kmap,
int i, int j, int k); int i, int j, int k);
PermeabilityKind fillTensor(Opm::EclipseStateConstPtr eclState,
std::vector<const std::vector<double>*>& tensor, typedef GridPropertyAccess::ArrayPolicy::ExtractFromDeck<double> PermArray;
std::array<int,9>& kmap);
struct PermTag {};
typedef GridPropertyAccess::Compressed<PermArray, PermTag> PermComponent;
PermComponent
extractPermComponent(EclipseStateConstPtr ecl,
const std::string& kw,
const int* global_cell);
PermeabilityKind
fillTensor(EclipseStateConstPtr eclState,
const int* global_cell,
std::vector<PermComponent>& tensor,
std::array<int,9>& kmap);
} // anonymous namespace } // anonymous namespace
@ -68,16 +86,15 @@ namespace Opm
void RockFromDeck::assignPorosity(Opm::EclipseStateConstPtr eclState, void RockFromDeck::assignPorosity(Opm::EclipseStateConstPtr eclState,
int number_of_cells, const int* global_cell) int number_of_cells, const int* global_cell)
{ {
porosity_.assign(number_of_cells, 1.0); typedef GridPropertyAccess::ArrayPolicy
if (eclState->hasDoubleGridProperty("PORO")) { ::ExtractFromDeck<double> Array;
const std::vector<double>& poro =
eclState->getDoubleGridProperty("PORO")->getData(); Array poro_glob(eclState, "PORO", 1.0);
for (int c = 0; c < int(porosity_.size()); ++c) { GridPropertyAccess::Compressed<Array> poro(poro_glob, global_cell);
const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
assert(0 <= c && c < (int) porosity_.size()); porosity_.clear(); porosity_.reserve(number_of_cells);
assert(0 <= deck_pos && deck_pos < (int) poro.size()); for (int c = 0; c < number_of_cells; ++c) {
porosity_[c] = poro[deck_pos]; porosity_.push_back(poro[c]);
}
} }
} }
@ -95,26 +112,23 @@ namespace Opm
permeability_.assign(dim * dim * nc, 0.0); permeability_.assign(dim * dim * nc, 0.0);
std::vector<const std::vector<double>*> tensor; std::vector<PermComponent> tensor;
tensor.reserve(10); tensor.reserve(6);
const std::vector<double> zero(num_global_cells, 0.0);
tensor.push_back(&zero);
std::array<int,9> kmap; std::array<int,9> kmap;
PermeabilityKind pkind = fillTensor(eclState, tensor, kmap); PermeabilityKind pkind = fillTensor(eclState, global_cell,
tensor, kmap);
if (pkind == Invalid) { if (pkind == Invalid) {
OPM_THROW(std::runtime_error, "Invalid permeability field."); OPM_THROW(std::runtime_error, "Invalid permeability field.");
} }
if (tensor.size() > 1) { assert (! tensor.empty());
const int* gc = global_cell; {
int off = 0; int off = 0;
for (int c = 0; c < nc; ++c, off += dim*dim) { for (int c = 0; c < nc; ++c, off += dim*dim) {
// SharedPermTensor K(dim, dim, &permeability_[off]); // SharedPermTensor K(dim, dim, &permeability_[off]);
int kix = 0; int kix = 0;
const int glob = (gc == NULL) ? c : gc[c];
for (int i = 0; i < dim; ++i) { for (int i = 0; i < dim; ++i) {
for (int j = 0; j < dim; ++j, ++kix) { for (int j = 0; j < dim; ++j, ++kix) {
@ -125,7 +139,8 @@ namespace Opm
// values in the resulting array are the same // values in the resulting array are the same
// in either order when viewed contiguously // in either order when viewed contiguously
// because fillTensor() enforces symmetry. // because fillTensor() enforces symmetry.
permeability_[off + (i + dim*j)] = (*tensor[kmap[kix]])[glob]; permeability_[off + (i + dim*j)] =
tensor[kmap[kix]][c];
} }
// K(i,i) = std::max(K(i,i), perm_threshold); // K(i,i) = std::max(K(i,i), perm_threshold);
@ -229,8 +244,8 @@ namespace Opm
void setScalarPermIfNeeded(std::array<int,9>& kmap, void setScalarPermIfNeeded(std::array<int,9>& kmap,
int i, int j, int k) int i, int j, int k)
{ {
if (kmap[j] == 0) { kmap[j] = kmap[i]; } if (kmap[j] < 0) { kmap[j] = kmap[i]; }
if (kmap[k] == 0) { kmap[k] = kmap[i]; } if (kmap[k] < 0) { kmap[k] = kmap[i]; }
} }
/// @brief /// @brief
@ -265,16 +280,20 @@ namespace Opm
/// ///
/// @param [out] tensor /// @param [out] tensor
/// @param [out] kmap /// @param [out] kmap
PermeabilityKind fillTensor(Opm::EclipseStateConstPtr eclState, PermeabilityKind
std::vector<const std::vector<double>*>& tensor, fillTensor(EclipseStateConstPtr eclState,
std::array<int,9>& kmap) const int* global_cell,
std::vector<PermComponent>& tensor,
std::array<int,9>& kmap)
{ {
PermeabilityKind kind = classifyPermeability(eclState); PermeabilityKind kind = classifyPermeability(eclState);
if (kind == Invalid) { if (kind == Invalid) {
OPM_THROW(std::runtime_error, "Invalid set of permeability fields given."); 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; } assert (tensor.empty());
for (int i = 0; i < 9; ++i) { kmap[i] = -1; }
enum { xx, xy, xz, // 0, 1, 2 enum { xx, xy, xz, // 0, 1, 2
yx, yy, yz, // 3, 4, 5 yx, yy, yz, // 3, 4, 5
@ -284,43 +303,53 @@ namespace Opm
// 1st row: [ kxx, kxy ], kxz handled in kzx // 1st row: [ kxx, kxy ], kxz handled in kzx
if (eclState->hasDoubleGridProperty("PERMX" )) { if (eclState->hasDoubleGridProperty("PERMX" )) {
kmap[xx] = tensor.size(); kmap[xx] = tensor.size();
tensor.push_back(&eclState->getDoubleGridProperty("PERMX")->getData()); tensor.push_back(extractPermComponent(eclState, "PERMX", global_cell));
setScalarPermIfNeeded(kmap, xx, yy, zz); setScalarPermIfNeeded(kmap, xx, yy, zz);
} }
if (eclState->hasDoubleGridProperty("PERMXY")) { {
kmap[xy] = kmap[yx] = tensor.size(); // Enforce symmetry. kmap[xy] = kmap[yx] = tensor.size(); // Enforce symmetry.
tensor.push_back(&eclState->getDoubleGridProperty("PERMXY")->getData()); tensor.push_back(extractPermComponent(eclState, "PERMXY", global_cell));
} }
// ----------------------------------------------------------- // -----------------------------------------------------------
// 2nd row: [ kyy, kyz ], kyx handled in kxy // 2nd row: [ kyy, kyz ], kyx handled in kxy
if (eclState->hasDoubleGridProperty("PERMY" )) { if (eclState->hasDoubleGridProperty("PERMY" )) {
kmap[yy] = tensor.size(); kmap[yy] = tensor.size();
tensor.push_back(&eclState->getDoubleGridProperty("PERMY")->getData()); tensor.push_back(extractPermComponent(eclState, "PERMY", global_cell));
setScalarPermIfNeeded(kmap, yy, zz, xx); setScalarPermIfNeeded(kmap, yy, zz, xx);
} }
if (eclState->hasDoubleGridProperty("PERMYZ")) { {
kmap[yz] = kmap[zy] = tensor.size(); // Enforce symmetry. kmap[yz] = kmap[zy] = tensor.size(); // Enforce symmetry.
tensor.push_back(&eclState->getDoubleGridProperty("PERMYZ")->getData()); tensor.push_back(extractPermComponent(eclState, "PERMYZ", global_cell));
} }
// ----------------------------------------------------------- // -----------------------------------------------------------
// 3rd row: [ kzx, kzz ], kzy handled in kyz // 3rd row: [ kzx, kzz ], kzy handled in kyz
if (eclState->hasDoubleGridProperty("PERMZX")) { {
kmap[zx] = kmap[xz] = tensor.size(); // Enforce symmetry. kmap[zx] = kmap[xz] = tensor.size(); // Enforce symmetry.
tensor.push_back(&eclState->getDoubleGridProperty("PERMZX")->getData()); tensor.push_back(extractPermComponent(eclState, "PERMZX", global_cell));
} }
if (eclState->hasDoubleGridProperty("PERMZ" )) { if (eclState->hasDoubleGridProperty("PERMZ" )) {
kmap[zz] = tensor.size(); kmap[zz] = tensor.size();
tensor.push_back(&eclState->getDoubleGridProperty("PERMZ")->getData()); tensor.push_back(extractPermComponent(eclState, "PERMZ", global_cell));
setScalarPermIfNeeded(kmap, zz, xx, yy); setScalarPermIfNeeded(kmap, zz, xx, yy);
} }
return kind; return kind;
} }
PermComponent
extractPermComponent(EclipseStateConstPtr ecl,
const std::string& kw,
const int* global_cell)
{
PermArray k(ecl, kw, 0.0); // return 0.0 if not present.
return PermComponent(k, global_cell);
}
} // anonymous namespace } // anonymous namespace
} // namespace Opm } // namespace Opm