Loading...
Searching...
No Matches
Rock_impl.hpp
/*
Copyright 2010 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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ROCK_IMPL_HEADER_INCLUDED
#define OPM_ROCK_IMPL_HEADER_INCLUDED
#include <fstream>
#include <array>
#include <opm/output/eclipse/EclipseGridInspector.hpp>
namespace Opm
{
// ----- Methods of Rock -----
template <int dim>
Rock<dim>::Rock()
: permeability_kind_(Invalid)
{
}
template <int dim>
void Rock<dim>::init(const Opm::Deck& deck,
const std::vector<int>& global_cell,
const double perm_threshold)
{
// This code is mostly copied from ReservoirPropertyCommon::init(...).
static_assert(dim == 3, "");
permfield_valid_.assign(global_cell.size(), false);
assignPorosity (deck, global_cell);
assignPermeability(deck, global_cell, perm_threshold);
}
template <int dim>
void Rock<dim>::init(const int num_cells,
const double uniform_poro,
const double uniform_perm)
{
permfield_valid_.assign(num_cells, true);
porosity_.assign(num_cells, uniform_poro);
permeability_.assign(dim*dim*num_cells, 0.0);
for (int i = 0; i < num_cells; ++i) {
SharedPermTensor K = permeabilityModifiable(i);
for (int dd = 0; dd < dim; ++dd) {
K(dd, dd) = uniform_perm;
}
}
}
template <int dim>
double Rock<dim>::porosity(int cell_index) const
{
return porosity_[cell_index];
}
template <int dim>
typename Rock<dim>::PermTensor
Rock<dim>::permeability(int cell_index) const
{
assert (permfield_valid_[cell_index]);
const PermTensor K(dim, dim, &permeability_[dim*dim*cell_index]);
return K;
}
template <int dim>
typename Rock<dim>::SharedPermTensor
Rock<dim>::permeabilityModifiable(int cell_index)
{
// Typically only used for assigning synthetic perm values.
SharedPermTensor K(dim, dim, &permeability_[dim*dim*cell_index]);
// Trust caller!
permfield_valid_[cell_index] = std::vector<unsigned char>::value_type(1);
return K;
}
// ------ Private methods ------
template <int dim>
void Rock<dim>::assignPorosity(const Opm::Deck& deck,
const std::vector<int>& global_cell)
{
porosity_.assign(global_cell.size(), 1.0);
if (deck.hasKeyword("PORO")) {
const std::vector<double>& poro = deck["PORO"].back().getSIDoubleData();
for (int c = 0; c < int(porosity_.size()); ++c) {
porosity_[c] = poro[global_cell[c]];
}
}
}
template <int dim>
void Rock<dim>::assignPermeability(const Opm::Deck& deck,
const std::vector<int>& global_cell,
double perm_threshold)
{
Opm::EclipseGridInspector insp(deck);
std::array<int, 3> dims = insp.gridSize();
int num_global_cells = dims[0]*dims[1]*dims[2];
assert (num_global_cells > 0);
permeability_.assign(dim * dim * global_cell.size(), 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);
static_assert(dim == 3, "");
std::array<int,9> kmap;
permeability_kind_ = fillTensor(deck, tensor, kmap);
// Assign permeability values only if such values are
// given in the input deck represented by 'deck'. 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 nc = global_cell.size();
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 = global_cell[c];
for (int i = 0; i < dim; ++i) {
for (int j = 0; j < dim; ++j, ++kix) {
K(i,j) = (*tensor[kmap[kix]])[glob];
}
K(i,i) = std::max(K(i,i), perm_threshold);
}
permfield_valid_[c] = std::vector<unsigned char>::value_type(1);
}
}
}
} // namespace Opm
#endif // OPM_ROCK_IMPL_HEADER_INCLUDED
Opm::Rock::permeability
PermTensor permeability(int cell_index) const
Read-access to permeability.
Definition Rock_impl.hpp:84
Opm::Rock::permeabilityModifiable
SharedPermTensor permeabilityModifiable(int cell_index)
Read- and write-access to permeability.
Definition Rock_impl.hpp:95
Opm::Rock::SharedPermTensor
SharedCMatrix SharedPermTensor
Tensor type for read and write access to permeability.
Definition Rock.hpp:43
Opm::Rock::PermTensor
ImmutableCMatrix PermTensor
Tensor type for read-only access to permeability.
Definition Rock.hpp:39
Opm::Rock::Rock
Rock()
Default constructor.
Definition Rock_impl.hpp:37
Opm::Rock::porosity
double porosity(int cell_index) const
Read-access to porosity.
Definition Rock_impl.hpp:76
Opm::Rock::init
void init(const Opm::Deck &deck, const std::vector< int > &global_cell, const double perm_threshold=0.0)
Initialize from a grdecl file.
Definition Rock_impl.hpp:44
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::zero
void zero(Matrix &A)
Zero-fill a.
Definition Matrix.hpp:602