OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

remove EclipseGridParser compatibility methods from all classes

Browse Source
This commit is contained in:
Andreas Lauser 2014-04-17 12:04:31 +02:00
parent abd341ab7f
commit f360562aee
24 changed files with 27 additions and 1718 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
opm/core/props/BlackoilPropertiesFromDeck.cpp
View File

@ -23,15 +23,6 @@
namespace Opm
{
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
bool init_rock)
{
init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids,
grid.dimensions, init_rock);
}
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
bool init_rock)
@ -40,15 +31,6 @@ namespace Opm
grid.cell_centroids, grid.dimensions, init_rock);
}
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,
bool init_rock)
{
init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids,
grid.dimensions, param, init_rock);
}
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,

82
opm/core/props/BlackoilPropertiesFromDeck.hpp
View File

@ -25,7 +25,6 @@
#include <opm/core/props/rock/RockFromDeck.hpp>
#include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
@ -42,14 +41,6 @@ namespace Opm
class BlackoilPropertiesFromDeck : public BlackoilPropertiesInterface
{
public:
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
/// \param[in] 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.
BlackoilPropertiesFromDeck(const EclipseGridParser &deck,
const UnstructuredGrid& grid, bool init_rock=true );
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
@ -58,22 +49,6 @@ namespace Opm
BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid, bool init_rock=true );
/// Initialize from deck, grid and parameters.
/// \param[in] deck Deck input parser
/// \param[in] 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.
/// \param[in] param Parameters. Accepted parameters include:
/// pvt_tab_size (200) number of uniform sample points for dead-oil pvt tables.
/// sat_tab_size (200) number of uniform sample points for saturation tables.
/// threephase_model("simple") three-phase relperm model (accepts "simple" and "stone2").
/// For both size parameters, a 0 or negative value indicates that no spline fitting is to
/// be done, and the input fluid data used directly for linear interpolation.
BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,
bool init_rock=true);
/// Initialize from deck, grid and parameters.
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
@ -90,41 +65,22 @@ namespace Opm
const parameter::ParameterGroup& param,
bool init_rock=true);
template<class T>
BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
bool init_rock=true);
template<class T>
BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock=true);
template<class T>
template<class CentroidIterator>
BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
const CentroidIterator& begin_cell_centroids,
int dimension,
bool init_rock=true);
template<class T>
template<class CentroidIterator>
BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
const CentroidIterator& begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock=true);
@ -250,38 +206,20 @@ namespace Opm
double* smax) const;
private:
template<class T>
void init(const EclipseGridParser& deck,
template<class CentroidIterator>
void init(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
const CentroidIterator& begin_cell_centroids,
int dimension,
bool init_rock);
template<class T>
void init(const EclipseGridParser& deck,
template<class CentroidIterator>
void init(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock);
template<class T>
void init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
bool init_rock);
template<class T>
void init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
const CentroidIterator& begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock);

12
opm/core/props/IncompPropertiesFromDeck.cpp
View File

@ -26,18 +26,6 @@
namespace Opm
{
IncompPropertiesFromDeck::IncompPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid)
{
rock_.init(deck, grid);
pvt_.init(deck);
satprops_.init(deck, grid, 200);
if (pvt_.numPhases() != satprops_.numPhases()) {
OPM_THROW(std::runtime_error, "IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
}
IncompPropertiesFromDeck::IncompPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid)
{

8
opm/core/props/IncompPropertiesFromDeck.hpp
View File

@ -45,14 +45,6 @@ namespace Opm
class IncompPropertiesFromDeck : public IncompPropertiesInterface
{
public:
/// Initialize from deck and grid.
/// \param deck Deck input parser
/// \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.
IncompPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid);
/// Initialize from deck and grid.
/// \param deck Deck input parser
/// \param grid Grid to which property object applies, needed for the

40
opm/core/props/phaseUsageFromDeck.hpp
View File

@ -20,9 +20,8 @@
#ifndef OPM_PHASEUSAGEFROMDECK_HEADER_INCLUDED
#define OPM_PHASEUSAGEFROMDECK_HEADER_INCLUDED
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@ -68,43 +67,6 @@ namespace Opm
return pu;
}
/// Looks at presence of WATER, OIL and GAS keywords in deck
/// to determine active phases.
inline PhaseUsage phaseUsageFromDeck(const EclipseGridParser& deck)
{
PhaseUsage pu;
std::fill(pu.phase_used, pu.phase_used + BlackoilPhases::MaxNumPhases, 0);
// Discover phase usage.
if (deck.hasField("WATER")) {
pu.phase_used[BlackoilPhases::Aqua] = 1;
}
if (deck.hasField("OIL")) {
pu.phase_used[BlackoilPhases::Liquid] = 1;
}
if (deck.hasField("GAS")) {
pu.phase_used[BlackoilPhases::Vapour] = 1;
}
pu.num_phases = 0;
for (int i = 0; i < BlackoilPhases::MaxNumPhases; ++i) {
pu.phase_pos[i] = pu.num_phases;
pu.num_phases += pu.phase_used[i];
}
// Only 2 or 3 phase systems handled.
if (pu.num_phases < 2 || pu.num_phases > 3) {
OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases.");
}
// We need oil systems, since we do not support the keywords needed for
// water-gas systems.
if (!pu.phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems.");
}
return pu;
}
/// Looks at presence of WATER, OIL and GAS keywords in deck
/// to determine active phases.
inline PhaseUsage phaseUsageFromDeck(Opm::DeckConstPtr newParserDeck)

55
opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
View File

@ -21,7 +21,6 @@
#include "config.h"
#include <opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
@ -34,60 +33,6 @@ namespace Opm
{
}
void PvtPropertiesIncompFromDeck::init(const EclipseGridParser& deck)
{
// If we need multiple regions, this class and the SinglePvt* classes must change.
int region_number = 0;
PhaseUsage phase_usage = phaseUsageFromDeck(deck);
if (phase_usage.phase_used[PhaseUsage::Vapour] ||
!phase_usage.phase_used[PhaseUsage::Aqua] ||
!phase_usage.phase_used[PhaseUsage::Liquid]) {
OPM_THROW(std::runtime_error, "PvtPropertiesIncompFromDeck::init() -- must have gas and oil phases (only) in deck input.\n");
}
// Surface densities. Accounting for different orders in eclipse and our code.
if (deck.hasField("DENSITY")) {
const std::vector<double>& d = deck.getDENSITY().densities_[region_number];
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
surface_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = d[ECL_water];
surface_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil];
} else {
OPM_THROW(std::runtime_error, "Input is missing DENSITY\n");
}
// Make reservoir densities the same as surface densities initially.
// We will modify them with formation volume factors if found.
reservoir_density_ = surface_density_;
// Water viscosity.
if (deck.hasField("PVTW")) {
const std::vector<double>& pvtw = deck.getPVTW().pvtw_[region_number];
if (pvtw[2] != 0.0 || pvtw[4] != 0.0) {
OPM_MESSAGE("Compressibility effects in PVTW are ignored.");
}
reservoir_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] /= pvtw[1];
viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3];
} else {
// Eclipse 100 default.
// viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = 0.5*Opm::prefix::centi*Opm::unit::Poise;
OPM_THROW(std::runtime_error, "Input is missing PVTW\n");
}
// Oil viscosity.
if (deck.hasField("PVCDO")) {
const std::vector<double>& pvcdo = deck.getPVCDO().pvcdo_[region_number];
if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) {
OPM_MESSAGE("Compressibility effects in PVCDO are ignored.");
}
reservoir_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] /= pvcdo[1];
viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3];
} else {
OPM_THROW(std::runtime_error, "Input is missing PVCDO\n");
}
}
void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr newParserDeck )
{
// If we need multiple regions, this class and the SinglePvt* classes must change.

4
opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp
View File

@ -20,7 +20,6 @@
#ifndef OPM_PVTPROPERTIESINCOMPFROMDECK_HEADER_INCLUDED
#define OPM_PVTPROPERTIESINCOMPFROMDECK_HEADER_INCLUDED
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <array>
@ -38,9 +37,6 @@ namespace Opm
/// Default constructor.
PvtPropertiesIncompFromDeck();
/// Initialize from deck.
void init(const EclipseGridParser& deck);
/// Initialize from deck.
void init(Opm::DeckConstPtr newParserDeck);

28
opm/core/props/rock/RockCompressibility.cpp
View File

@ -19,7 +19,6 @@
#include "config.h"
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
@ -41,33 +40,6 @@ namespace Opm
rock_comp_ = param.getDefault("rock_compressibility", 0.0)/unit::barsa;
}
RockCompressibility::RockCompressibility(const EclipseGridParser& deck)
: pref_(0.0),
rock_comp_(0.0)
{
if (deck.hasField("ROCKTAB")) {
const table_t& rt = deck.getROCKTAB().rocktab_;
if (rt.size() != 1) {
OPM_THROW(std::runtime_error, "Can only handle a single region in ROCKTAB.");
}
const int n = rt[0][0].size();
p_.resize(n);
poromult_.resize(n);
transmult_.resize(n);
for (int i = 0; i < n; ++i) {
p_[i] = rt[0][0][i];
poromult_[i] = rt[0][1][i];
transmult_[i] = rt[0][2][i];
}
} else if (deck.hasField("ROCK")) {
const ROCK& r = deck.getROCK();
pref_ = r.rock_compressibilities_[0][0];
rock_comp_ = r.rock_compressibilities_[0][1];
} else {
std::cout << "**** warning: no rock compressibility data found in deck (ROCK or ROCKTAB)." << std::endl;
}
}
RockCompressibility::RockCompressibility(Opm::DeckConstPtr newParserDeck)
: pref_(0.0),
rock_comp_(0.0)

5
opm/core/props/rock/RockCompressibility.hpp
View File

@ -27,16 +27,11 @@
namespace Opm
{
class EclipseGridParser;
namespace parameter { class ParameterGroup; }
class RockCompressibility
{
public:
/// Construct from input deck.
/// 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);

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

@ -21,6 +21,7 @@
#include "config.h"
#include <opm/core/props/rock/RockFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
@ -34,12 +35,8 @@ namespace Opm
{
enum PermeabilityKind { ScalarPerm, DiagonalPerm, TensorPerm, None, Invalid };
PermeabilityKind classifyPermeability(const EclipseGridParser& parser);
void setScalarPermIfNeeded(std::array<int,9>& kmap,
int i, int j, int k);
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);
@ -56,34 +53,6 @@ namespace Opm
{
}
/// Initialize from deck and cell mapping.
/// \param deck Deck input parser
/// \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 RockFromDeck::init(const EclipseGridParser& deck,
const UnstructuredGrid& grid)
{
init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims);
}
/// Initialize from deck and cell mapping.
/// \param deck Deck input parser
/// \param number_of_cells The number of cells in the grid.
/// \param global_cell The mapping fom local to global cell indices.
/// global_cell[i] is the corresponding global index of i.
/// \param cart_dims The size of the underlying cartesian grid.
void RockFromDeck::init(const EclipseGridParser& deck,
int number_of_cells, const int* global_cell,
const int* cart_dims)
{
assignPorosity(deck, number_of_cells, global_cell);
permfield_valid_.assign(number_of_cells, false);
const double perm_threshold = 0.0; // Maybe turn into parameter?
assignPermeability(deck, number_of_cells, global_cell, cart_dims, perm_threshold);
}
void RockFromDeck::init(Opm::DeckConstPtr newParserDeck,
int number_of_cells, const int* global_cell,
const int* cart_dims)
@ -95,20 +64,6 @@ namespace Opm
perm_threshold);
}
void RockFromDeck::assignPorosity(const EclipseGridParser& parser,
int number_of_cells, const int* global_cell)
{
porosity_.assign(number_of_cells, 1.0);
if (parser.hasField("PORO")) {
const std::vector<double>& poro = parser.getFloatingPointValue("PORO");
for (int c = 0; c < int(porosity_.size()); ++c) {
const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
porosity_[c] = poro[deck_pos];
}
}
}
void RockFromDeck::assignPorosity(Opm::DeckConstPtr newParserDeck,
int number_of_cells, const int* global_cell)
{
@ -124,61 +79,6 @@ namespace Opm
}
}
void RockFromDeck::assignPermeability(const EclipseGridParser& parser,
int number_of_cells,
const int* global_cell,
const int* cartdims,
double perm_threshold)
{
const int dim = 3;
const int num_global_cells = cartdims[0]*cartdims[1]*cartdims[2];
assert(num_global_cells > 0);
permeability_.assign(dim * dim * number_of_cells, 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(parser, 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 'parser'. 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) {
int off = 0;
for (int c = 0; c < number_of_cells; ++c, off += dim*dim) {
// SharedPermTensor K(dim, dim, &permeability_[off]);
int kix = 0;
const int glob = (global_cell == NULL) ? c : 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];
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);
}
}
}
void RockFromDeck::assignPermeability(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
@ -236,73 +136,6 @@ namespace Opm
}
namespace {
/// @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 parser [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(const EclipseGridParser& parser)
{
const bool xx = parser.hasField("PERMX" );
const bool xy = parser.hasField("PERMXY");
const bool xz = parser.hasField("PERMXZ");
const bool yx = parser.hasField("PERMYX");
const bool yy = parser.hasField("PERMY" );
const bool yz = parser.hasField("PERMYZ");
const bool zx = parser.hasField("PERMZX");
const bool zy = parser.hasField("PERMZY");
const bool zz = parser.hasField("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
/// Classify and verify a given permeability specification
/// from a structural point of view. In particular, we
@ -397,107 +230,6 @@ namespace Opm
if (kmap[k] == 0) { kmap[k] = kmap[i]; }
}
/// @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(const EclipseGridParser& parser,
std::vector<const std::vector<double>*>& tensor,
std::array<int,9>& kmap)
{
PermeabilityKind kind = classifyPermeability(parser);
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 (parser.hasField("PERMX" )) {
kmap[xx] = tensor.size();
tensor.push_back(&parser.getFloatingPointValue("PERMX" ));
setScalarPermIfNeeded(kmap, xx, yy, zz);
}
if (parser.hasField("PERMXY")) {
kmap[xy] = kmap[yx] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMXY"));
}
if (parser.hasField("PERMXZ")) {
kmap[xz] = kmap[zx] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMXZ"));
}
// -----------------------------------------------------------
// 2nd row: [kyx, kyy, kyz]
if (parser.hasField("PERMYX")) {
kmap[yx] = kmap[xy] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMYX"));
}
if (parser.hasField("PERMY" )) {
kmap[yy] = tensor.size();
tensor.push_back(&parser.getFloatingPointValue("PERMY" ));
setScalarPermIfNeeded(kmap, yy, zz, xx);
}
if (parser.hasField("PERMYZ")) {
kmap[yz] = kmap[zy] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMYZ"));
}
// -----------------------------------------------------------
// 3rd row: [kzx, kzy, kzz]
if (parser.hasField("PERMZX")) {
kmap[zx] = kmap[xz] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMZX"));
}
if (parser.hasField("PERMZY")) {
kmap[zy] = kmap[yz] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMZY"));
}
if (parser.hasField("PERMZ" )) {
kmap[zz] = tensor.size();
tensor.push_back(&parser.getFloatingPointValue("PERMZ" ));
setScalarPermIfNeeded(kmap, zz, xx, yy);
}
return kind;
}
/// @brief
/// Extract pointers to appropriate tensor components from
/// input deck. The permeability tensor is, generally,

28
opm/core/props/rock/RockFromDeck.hpp
View File

@ -20,9 +20,6 @@
#ifndef OPM_ROCKFROMDECK_HEADER_INCLUDED
#define OPM_ROCKFROMDECK_HEADER_INCLUDED
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <vector>
@ -38,24 +35,7 @@ namespace Opm
/// Default constructor.
RockFromDeck();
/// Initialize from deck and grid.
/// \param deck Deck input parser
/// \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(const EclipseGridParser& deck,
const UnstructuredGrid& grid);
/// Initialize from deck and cell mapping.
/// \param deck Deck input parser
/// \param number_of_cells The number of cells in the grid.
/// \param global_cell The mapping fom local to global cell indices.
/// global_cell[i] is the corresponding global index of i.
/// \param cart_dims The size of the underlying cartesian grid.
void init(const EclipseGridParser& deck,
int number_of_cells, const int* global_cell,
const int* cart_dims);
/// Initialize from deck and cell mapping.
/// \param newParserDeck Deck produced by the opm-parser code
/// \param number_of_cells The number of cells in the grid.
/// \param global_cell The mapping fom local to global cell indices.
@ -92,17 +72,9 @@ namespace Opm
}
private:
void assignPorosity(const EclipseGridParser& parser,
int number_of_cells,
const int* global_cell);
void assignPorosity(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell);
void assignPermeability(const EclipseGridParser& parser,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
const double perm_threshold);
void assignPermeability(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,

53
opm/core/props/satfunc/SaturationPropsFromDeck.hpp
View File

@ -22,7 +22,6 @@
#include <opm/core/props/satfunc/SaturationPropsInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SatFuncStone2.hpp>
#include <opm/core/props/satfunc/SatFuncSimple.hpp>
@ -52,17 +51,6 @@ namespace Opm
/// Default constructor.
SaturationPropsFromDeck();
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
/// \param[in] 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.
/// \param[in] samples Number of uniform sample points for saturation tables.
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
void init(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const int samples);
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
@ -74,29 +62,11 @@ namespace Opm
const UnstructuredGrid& grid,
const int samples);
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
/// \param[in] newParserDeck Deck input parser
/// \param[in] number_of_cells The number of cells of the 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.
/// \param[in] global_cell The mapping from local cell indices of the grid to
/// global cell indices used in the deck.
/// \param[in] begin_cell_centroids Pointer to the first cell_centroid of the grid.
/// \param[in] dimensions The dimensions of the grid.
/// \param[in] samples Number of uniform sample points for saturation tables.
/// \tparam T The iterator Type for the cell centroids.
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
template<class T>
void init(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const int samples);
/// grid) to logical cartesian indices consistent with the
/// deck.
/// \param[in] global_cell The mapping from local cell indices of the grid to
@ -178,27 +148,6 @@ namespace Opm
typedef SatFuncSet Funcs;
const Funcs& funcForCell(const int cell) const;
template<class T>
void initEPS(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSHyst(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSKey(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam);
template<class T>
void initEPS(Opm::DeckConstPtr newParserDeck,
int number_of_cells,

463
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
View File

@ -47,126 +47,15 @@ namespace Opm
{
}
/// Initialize from deck.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::init(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const int samples)
{
init(deck, grid.number_of_cells, grid.global_cell, grid.cell_centroids,
grid.dimensions, samples);
}
/// Initialize from deck.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
const int samples)
{
init(newParserDeck, grid.number_of_cells,
grid.global_cell, grid.cell_centroids,
grid.dimensions, samples);
}
/// Initialize from deck.
template <class SatFuncSet>
template< class T>
void SaturationPropsFromDeck<SatFuncSet>::init(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const int samples)
{
phase_usage_ = phaseUsageFromDeck(deck);
// Extract input data.
// Oil phase should be active.
if (!phase_usage_.phase_used[Liquid]) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- oil phase must be active.");
}
// Obtain SATNUM, if it exists, and create cell_to_func_.
// Otherwise, let the cell_to_func_ mapping be just empty.
int satfuncs_expected = 1;
if (deck.hasField("SATNUM")) {
const std::vector<int>& satnum = deck.getIntegerValue("SATNUM");
satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
cell_to_func_.resize(number_of_cells);
for (int cell = 0; cell < number_of_cells; ++cell) {
const int deck_pos = (global_cell == NULL) ? cell : global_cell[cell];
cell_to_func_[cell] = satnum[deck_pos] - 1;
}
}
// Find number of tables, check for consistency.
enum { Uninitialized = -1 };
int num_tables = Uninitialized;
if (phase_usage_.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
num_tables = swof_table.size();
if (num_tables < satfuncs_expected) {
OPM_THROW(std::runtime_error, "Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected);
}
}
if (phase_usage_.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
int num_sgof_tables = sgof_table.size();
if (num_sgof_tables < satfuncs_expected) {
OPM_THROW(std::runtime_error, "Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected);
}
if (num_tables == Uninitialized) {
num_tables = num_sgof_tables;
} else if (num_tables != num_sgof_tables) {
OPM_THROW(std::runtime_error, "Inconsistent number of tables in SWOF and SGOF.");
}
}
// Initialize tables.
satfuncset_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(deck, table, phase_usage_, samples);
}
// Saturation table scaling
do_hyst_ = false;
do_eps_ = false;
do_3pt_ = false;
if (deck.hasField("ENDSCALE")) {
//if (!phase_usage_.phase_used[Aqua] || !phase_usage_.phase_used[Liquid] || phase_usage_.phase_used[Vapour]) {
// OPM_THROW(std::runtime_error, "Currently endpoint-scaling limited to oil-water systems without gas.");
//}
if (deck.getENDSCALE().dir_switch_ != std::string("NODIR")) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'NODIR' accepted.");
}
if (deck.getENDSCALE().revers_switch_ != std::string("REVERS")) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'REVERS' accepted.");
}
if (deck.hasField("SCALECRS")) {
if (deck.getSCALECRS().scalecrs_ == std::string("YES")) {
do_3pt_ = true;
}
}
do_eps_ = true;
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions);
// For now, a primitive detection of hysteresis. TODO: SATOPTS HYSTER/ and EHYSTR
do_hyst_ = deck.hasField("ISWL") || deck.hasField("ISWU") || deck.hasField("ISWCR") || deck.hasField("ISGL") ||
deck.hasField("ISGU") || deck.hasField("ISGCR") || deck.hasField("ISOWCR") || deck.hasField("ISOGCR");
if (do_hyst_) {
if (deck.hasField("KRW") || deck.hasField("KRG") || deck.hasField("KRO") || deck.hasField("KRWR") ||
deck.hasField("KRGR") || deck.hasField("KRORW") || deck.hasField("KRORG") ||
deck.hasField("IKRW") || deck.hasField("IKRG") || deck.hasField("IKRO") || deck.hasField("IKRWR") ||
deck.hasField("IKRGR") || deck.hasField("IKRORW") || deck.hasField("IKRORG") ) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Currently hysteresis and relperm value scaling can not be combined.");
}
initEPSHyst(deck, number_of_cells, global_cell, begin_cell_centroids,
dimensions);
}
//OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Under construction ...");
}
this->init(newParserDeck, grid.number_of_cells,
grid.global_cell, grid.cell_centroids,
grid.dimensions, samples);
}
/// Initialize from deck.
@ -532,89 +421,6 @@ namespace Opm
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
// Initialize saturation scaling parameters
template <class SatFuncSet>
template <class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
int dimensions)
{
std::vector<double> swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr;
std::vector<double> krw, krg, kro, krwr, krgr, krorw, krorg;
// Initialize saturation scaling parameter
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWL"), swl);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWU"), swu);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWCR"), swcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGL"), sgl);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGU"), sgu);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGCR"), sgcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SOWCR"), sowcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SOGCR"), sogcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRW"), krw);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRG"), krg);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRO"), kro);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRWR"), krwr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRGR"), krgr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRORW"), krorw);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRORG"), krorg);
eps_transf_.resize(number_of_cells);
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour];
const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
for (int cell = 0; cell < number_of_cells; ++cell) {
if (oilWater) {
// ### krw
initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos],
funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw);
// ### krow
initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro);
} else if (oilGas) {
// ### krg
initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos],
funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg);
// ### krog
initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro);
} else if (threephase) {
// ### krw
initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_,
funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw);
// ### krow
initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_,
funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro);
// ### krg
initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_,
funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg);
// ### krog
initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_,
funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro);
}
}
}
// Initialize saturation scaling parameters
template <class SatFuncSet>
template<class T>
@ -698,91 +504,6 @@ namespace Opm
}
}
// Initialize hysteresis saturation scaling parameters
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
int dimensions)
{
std::vector<double> iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr;
std::vector<double> ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg;
// Initialize hysteresis saturation scaling parameters
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWL"), iswl);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWU"), iswu);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWCR"), iswcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGL"), isgl);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGU"), isgu);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGCR"), isgcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISOWCR"), isowcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISOGCR"), isogcr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRW"), ikrw);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRG"), ikrg);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRO"), ikro);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRWR"), ikrwr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRGR"), ikrgr);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRORW"), ikrorw);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRORG"), ikrorg);
eps_transf_hyst_.resize(number_of_cells);
sat_hyst_.resize(number_of_cells);
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour];
const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
for (int cell = 0; cell < number_of_cells; ++cell) {
if (oilWater) {
// ### krw
initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos],
funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw);
// ### krow
initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro);
} else if (oilGas) {
// ### krg
initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos],
funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg);
// ### krog
initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro);
} else if (threephase) {
// ### krw
initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_,
funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw);
// ### krow
initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_,
funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro);
// ### krg
initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_,
funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg);
// ### krog
initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_,
funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro);
}
}
}
// Initialize hysteresis saturation scaling parameters
template <class SatFuncSet>
template<class T>
@ -867,184 +588,6 @@ namespace Opm
}
}
// Initialize saturation scaling parameter
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam)
{
const bool useAqua = phase_usage_.phase_used[Aqua];
const bool useLiquid = phase_usage_.phase_used[Liquid];
const bool useVapour = phase_usage_.phase_used[Vapour];
bool useKeyword = deck.hasField(keyword);
bool hasENPTVD = deck.hasField("ENPTVD");
bool hasENKRVD = deck.hasField("ENKRVD");
int itab = 0;
std::vector<std::vector<std::vector<double> > > table_dummy;
std::vector<std::vector<std::vector<double> > >& table = table_dummy;
// Active keyword assigned default values for each cell (in case of possible box-wise assignment)
int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua];
int phase_pos_vapour = phase_usage_.phase_pos[BlackoilPhases::Vapour];
if ((keyword[0] == 'S' && (useKeyword || hasENPTVD)) || (keyword[1] == 'S' && useKeyword) ) {
if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) {
if (useAqua && (useKeyword || deck.getENPTVD().mask_[0])) {
itab = 1;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smin_[phase_pos_aqua];
}
} else if (keyword == std::string("SWCR") || keyword == std::string("ISWCR") ) {
if (useAqua && (useKeyword || deck.getENPTVD().mask_[1])) {
itab = 2;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).swcr_;
}
} else if (keyword == std::string("SWU") || keyword == std::string("ISWU") ) {
if (useAqua && (useKeyword || deck.getENPTVD().mask_[2])) {
itab = 3;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smax_[phase_pos_aqua];
}
} else if (keyword == std::string("SGL") || keyword == std::string("ISGL") ) {
if (useVapour && (useKeyword || deck.getENPTVD().mask_[3])) {
itab = 4;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smin_[phase_pos_vapour];
}
} else if (keyword == std::string("SGCR") || keyword == std::string("ISGCR") ) {
if (useVapour && (useKeyword || deck.getENPTVD().mask_[4])) {
itab = 5;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sgcr_;
}
} else if (keyword == std::string("SGU") || keyword == std::string("ISGU") ) {
if (useVapour && (useKeyword || deck.getENPTVD().mask_[5])) {
itab = 6;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smax_[phase_pos_vapour];
}
} else if (keyword == std::string("SOWCR") || keyword == std::string("ISOWCR") ) {
if (useAqua && (useKeyword || deck.getENPTVD().mask_[6])) {
itab = 7;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sowcr_;
}
} else if (keyword == std::string("SOGCR") || keyword == std::string("ISOGCR") ) {
if (useVapour && (useKeyword || deck.getENPTVD().mask_[7])) {
itab = 8;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sogcr_;
}
} else {
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
table = deck.getENPTVD().table_;
}
} else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) {
if (keyword == std::string("KRW") || keyword == std::string("IKRW") ) {
if (useAqua && (useKeyword || deck.getENKRVD().mask_[0])) {
itab = 1;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwmax_;
}
} else if (keyword == std::string("KRG") || keyword == std::string("IKRG") ) {
if (useVapour && (useKeyword || deck.getENKRVD().mask_[1])) {
itab = 2;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgmax_;
}
} else if (keyword == std::string("KRO") || keyword == std::string("IKRO") ) {
if (useLiquid && (useKeyword || deck.getENKRVD().mask_[2])) {
itab = 3;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).kromax_;
}
} else if (keyword == std::string("KRWR") || keyword == std::string("IKRWR") ) {
if (useAqua && (useKeyword || deck.getENKRVD().mask_[3])) {
itab = 4;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwr_;
}
} else if (keyword == std::string("KRGR") || keyword == std::string("IKRGR") ) {
if (useVapour && (useKeyword || deck.getENKRVD().mask_[4])) {
itab = 5;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgr_;
}
} else if (keyword == std::string("KRORW") || keyword == std::string("IKRORW") ) {
if (useAqua && (useKeyword || deck.getENKRVD().mask_[5])) {
itab = 6;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorw_;
}
} else if (keyword == std::string("KRORG") || keyword == std::string("IKRORG") ) {
if (useVapour && (useKeyword || deck.getENKRVD().mask_[6])) {
itab = 7;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorg_;
}
} else {
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
table = deck.getENKRVD().table_;
}
}
if (scaleparam.empty()) {
return;
} else if (useKeyword) {
// Keyword values from deck
std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
const std::vector<double>& val = deck.getFloatingPointValue(keyword);
for (int c = 0; c < int(scaleparam.size()); ++c) {
const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
scaleparam[c] = val[deck_pos];
}
} else {
std::cout << "--- Scaling parameter '" << keyword << "' assigned via ";
if (keyword[0] == 'S')
deck.getENPTVD().write(std::cout);
else
deck.getENKRVD().write(std::cout);
for (int cell = 0; cell < number_of_cells; ++cell) {
int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell];
if (table[itab][jtab][0] != -1.0) {
std::vector<double>& depth = table[0][jtab];
std::vector<double>& val = table[itab][jtab];
double zc = UgGridHelpers
::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell,
dimensions),
dimensions-1);
if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval
scaleparam[cell] = linearInterpolation(depth, val, zc);
}
}
}
}
}
// Initialize saturation scaling parameter
template <class SatFuncSet>
template<class T>

49
opm/core/simulator/initState.hpp
View File

@ -20,7 +20,7 @@
#ifndef OPM_INITSTATE_HEADER_INCLUDED
#define OPM_INITSTATE_HEADER_INCLUDED
#include <opm/parser/eclipse/Deck/Deck.hpp> // DeckConstPtr
#include <opm/parser/eclipse/Deck/Deck.hpp>
struct UnstructuredGrid;
@ -28,7 +28,6 @@ namespace Opm
{
namespace parameter { class ParameterGroup; }
class EclipseGridParser;
class IncompPropertiesInterface;
class BlackoilPropertiesInterface;
@ -168,27 +167,7 @@ namespace Opm
template <class Props, class State>
void initStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
const EclipseGridParser& deck,
const double gravity,
State& state);
/// Initialize a two-phase state from input deck.
/// If EQUIL is present:
/// - saturation is set according to the water-oil contact,
/// - pressure is set to hydrostatic equilibrium.
/// Otherwise:
/// - saturation is set according to SWAT,
/// - pressure is set according to PRESSURE.
template <class FaceCells, class FCI, class CCI, class Props, class State>
void initStateFromDeck(int number_of_cells,
const int* global_cell,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const Props& props,
const EclipseGridParser& deck,
Opm::DeckConstPtr newParserDeck,
const double gravity,
State& state);
@ -203,31 +182,9 @@ namespace Opm
template <class Props, class State>
void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
const EclipseGridParser& deck,
Opm::DeckConstPtr newParserDeck,
const double gravity,
State& state);
/// Initialize a two-phase water-oil blackoil state from input deck.
/// If EQUIL is present:
/// - saturation is set according to the water-oil contact,
/// - pressure is set to hydrostatic equilibrium.
/// Otherwise:
/// - saturation is set according to SWAT,
/// - pressure is set according to PRESSURE.
/// In addition, this function sets surfacevol.
template <class FaceCells, class FCI, class CCI, class Props, class State>
void initBlackoilStateFromDeck(int number_of_cells,
const int* global_cell,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const Props& props,
const EclipseGridParser& deck,
const double gravity,
State& state);
/// Initialize a blackoil state from input deck.
template <class FaceCells, class FCI, class CCI, class Props, class State>
void initBlackoilStateFromDeck(int number_of_cells,

254
opm/core/simulator/initStateEquil.hpp
View File

@ -22,7 +22,6 @@
#include <opm/core/simulator/EquilibrationHelpers.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/RegionMapping.hpp>
@ -60,13 +59,6 @@ namespace Opm
* \param[in] deck Simulation deck, used to obtain EQUIL and related data.
* \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
*/
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const EclipseGridParser& deck,
const double gravity,
BlackoilState& state);
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr newParserDeck,
@ -188,54 +180,6 @@ namespace Opm
const std::vector<double> gas_saturation);
namespace DeckDependent {
inline
std::vector<EquilRecord>
getEquil(const EclipseGridParser& deck)
{
if (deck.hasField("EQUIL")) {
const EQUIL& eql = deck.getEQUIL();
typedef std::vector<EquilLine>::size_type sz_t;
const sz_t nrec = eql.equil.size();
std::vector<EquilRecord> ret;
ret.reserve(nrec);
for (sz_t r = 0; r < nrec; ++r) {
const EquilLine& rec = eql.equil[r];
EquilRecord record =
{
{ rec.datum_depth_ ,
rec.datum_depth_pressure_ }
,
{ rec.water_oil_contact_depth_ ,
rec.oil_water_cap_pressure_ }
,
{ rec.gas_oil_contact_depth_ ,
rec.gas_oil_cap_pressure_ }
,
rec.live_oil_table_index_
,
rec.wet_gas_table_index_
,
rec.N_
};
if (record.N != 0) {
OPM_THROW(std::domain_error,
"kw EQUIL, item 9: Only N=0 supported.");
}
ret.push_back(record);
}
return ret;
}
else {
OPM_THROW(std::domain_error,
"Deck does not provide equilibration data.");
}
}
inline
std::vector<EquilRecord>
getEquil(const Opm::DeckConstPtr newParserDeck)
@ -282,29 +226,6 @@ namespace Opm
}
}
inline
std::vector<int>
equilnum(const EclipseGridParser& deck,
const UnstructuredGrid& G )
{
std::vector<int> eqlnum;
if (deck.hasField("EQLNUM")) {
const std::vector<int>& e = deck.getIntegerValue("EQLNUM");
eqlnum.reserve(e.size());
std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
std::bind2nd(std::minus<int>(), 1));
}
else {
// No explicit equilibration region.
// All cells in region zero.
eqlnum.assign(G.number_of_cells, 0);
}
return eqlnum;
}
inline
std::vector<int>
equilnum(const Opm::DeckConstPtr newParserDeck,
@ -331,180 +252,7 @@ namespace Opm
}
template <class InputDeck>
class InitialStateComputer;
template <>
class InitialStateComputer<Opm::EclipseGridParser> {
public:
InitialStateComputer(const BlackoilPropertiesInterface& props,
const EclipseGridParser& deck ,
const UnstructuredGrid& G ,
const double grav = unit::gravity)
: pp_(props.numPhases(),
std::vector<double>(G.number_of_cells)),
sat_(props.numPhases(),
std::vector<double>(G.number_of_cells)),
rs_(G.number_of_cells),
rv_(G.number_of_cells)
{
// Get the equilibration records.
const std::vector<EquilRecord> rec = getEquil(deck);
// Create (inverse) region mapping.
const RegionMapping<> eqlmap(equilnum(deck, G));
// Create Rs functions.
rs_func_.reserve(rec.size());
if (deck.hasField("DISGAS")) {
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].live_oil_table_index > 0) {
if (deck.hasField("RSVD")) {
// TODO When this kw is actually parsed, also check for proper number of available tables
// For now, just use dummy ...
std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
std::vector<double> rs; rs.push_back(0.0); rs.push_back(100.0);
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, depth, rs));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
}
} else {
if (rec[i].goc.depth != rec[i].main.depth) {
OPM_THROW(std::runtime_error,
"Cannot initialise: when no explicit RSVD table is given, \n"
"datum depth must be at the gas-oil-contact. "
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
const double p_contact = rec[i].main.press;
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact));
}
}
} else {
for (size_t i = 0; i < rec.size(); ++i) {
rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
}
}
rv_func_.reserve(rec.size());
if (deck.hasField("VAPOIL")) {
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].wet_gas_table_index > 0) {
if (deck.hasField("RVVD")) {
// TODO When this kw is actually parsed, also check for proper number of available tables
// For now, just use dummy ...
std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
std::vector<double> rv; rv.push_back(0.0); rv.push_back(0.0001);
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props, cell, depth, rv));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
}
} else {
if (rec[i].goc.depth != rec[i].main.depth) {
OPM_THROW(std::runtime_error,
"Cannot initialise: when no explicit RVVD table is given, \n"
"datum depth must be at the gas-oil-contact. "
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
const double p_contact = rec[i].main.press + rec[i].goc.press;
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact));
}
}
} else {
for (size_t i = 0; i < rec.size(); ++i) {
rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
}
}
// Compute pressures, saturations, rs and rv factors.
calcPressSatRsRv(eqlmap, rec, props, G, grav);
// Modify oil pressure in no-oil regions so that the pressures of present phases can
// be recovered from the oil pressure and capillary relations.
}
typedef std::vector<double> Vec;
typedef std::vector<Vec> PVec; // One per phase.
const PVec& press() const { return pp_; }
const PVec& saturation() const { return sat_; }
const Vec& rs() const { return rs_; }
const Vec& rv() const { return rv_; }
private:
typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
typedef EquilReg<RhoCalc> EqReg;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
PVec pp_;
PVec sat_;
Vec rs_;
Vec rv_;
template <class RMap>
void
calcPressSatRsRv(const RMap& reg ,
const std::vector< EquilRecord >& rec ,
const Opm::BlackoilPropertiesInterface& props,
const UnstructuredGrid& G ,
const double grav)
{
typedef Miscibility::NoMixing NoMix;
for (typename RMap::RegionId
r = 0, nr = reg.numRegions();
r < nr; ++r)
{
const typename RMap::CellRange cells = reg.cells(r);
const int repcell = *cells.begin();
const RhoCalc calc(props, repcell);
const EqReg eqreg(rec[r], calc,
rs_func_[r], rv_func_[r],
props.phaseUsage());
PVec press = phasePressures(G, eqreg, cells, grav);
const PVec sat = phaseSaturations(eqreg, cells, props, press);
const int np = props.numPhases();
for (int p = 0; p < np; ++p) {
copyFromRegion(press[p], cells, pp_[p]);
copyFromRegion(sat[p], cells, sat_[p]);
}
if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
copyFromRegion(rs, cells, rs_);
copyFromRegion(rv, cells, rv_);
}
}
}
template <class CellRangeType>
void copyFromRegion(const Vec& source,
const CellRangeType& cells,
Vec& destination)
{
auto s = source.begin();
auto c = cells.begin();
const auto e = cells.end();
for (; c != e; ++c, ++s) {
destination[*c] = *s;
}
}
};
template <>
class InitialStateComputer<Opm::DeckConstPtr> {
class InitialStateComputer {
public:
InitialStateComputer(const BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr newParserDeck,

22
opm/core/simulator/initStateEquil_impl.hpp
View File

@ -735,33 +735,13 @@ namespace Opm
* \param[in] deck Simulation deck, used to obtain EQUIL and related data.
* \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
*/
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const EclipseGridParser& deck,
const double gravity,
BlackoilState& state)
{
typedef Equil::DeckDependent::InitialStateComputer<EclipseGridParser> ISC;
ISC isc(props, deck, grid, gravity);
const auto pu = props.phaseUsage();
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
? pu.phase_pos[BlackoilPhases::Liquid]
: pu.phase_pos[BlackoilPhases::Aqua];
state.pressure() = isc.press()[ref_phase];
state.saturation() = convertSats(isc.saturation());
state.gasoilratio() = isc.rs();
state.rv() = isc.rv();
// TODO: state.surfacevol() must be computed from s, rs, rv.
}
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr newParserDeck,
const double gravity,
BlackoilState& state)
{
typedef Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> ISC;
typedef Equil::DeckDependent::InitialStateComputer ISC;
ISC isc(props, newParserDeck, grid, gravity);
const auto pu = props.phaseUsage();
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]

179
opm/core/simulator/initState_impl.hpp
View File

@ -24,7 +24,6 @@
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/MonotCubicInterpolator.hpp>
#include <opm/core/utility/Units.hpp>
@ -673,129 +672,6 @@ namespace Opm
begin_cell_centroids, state);
}
/// Initialize a state from input deck.
template <class Props, class State>
void initStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
const EclipseGridParser& deck,
const double gravity,
State& state)
{
initStateFromDeck(grid.number_of_cells, grid.global_cell,
grid.number_of_faces, UgGridHelpers::faceCells(grid),
grid.face_centroids, grid.cell_centroids,
grid.dimensions, props, deck,
gravity, state);
}
template <class FaceCells, class FCI, class CCI, class Props, class State>
void initStateFromDeck(int number_of_cells,
const int* global_cell,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const Props& props,
const EclipseGridParser& deck,
const double gravity,
State& state)
{
const int num_phases = props.numPhases();
const PhaseUsage pu = phaseUsageFromDeck(deck);
if (num_phases != pu.num_phases) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): user specified property object with " << num_phases << " phases, "
"found " << pu.num_phases << " phases in deck.");
}
if (deck.hasField("EQUIL") && deck.hasField("PRESSURE")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
"condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
}
state.init(number_of_cells, number_of_faces, num_phases);
if (deck.hasField("EQUIL")) {
if (num_phases != 2) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas).");
}
// Set saturations depending on oil-water contact.
const EQUIL& equil= deck.getEQUIL();
if (equil.equil.size() != 1) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet.");
}
const double woc = equil.equil[0].water_oil_contact_depth_;
initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
props, woc, WaterBelow, state);
// Set pressure depending on densities and depths.
const double datum_z = equil.equil[0].datum_depth_;
const double datum_p = equil.equil[0].datum_depth_pressure_;
initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
props, woc, gravity, datum_z, datum_p, state);
} else if (deck.hasField("PRESSURE")) {
// Set saturations from SWAT/SGAS, pressure from PRESSURE.
std::vector<double>& s = state.saturation();
std::vector<double>& p = state.pressure();
const std::vector<double>& p_deck = deck.getFloatingPointValue("PRESSURE");
if (num_phases == 2) {
if (!pu.phase_used[BlackoilPhases::Aqua]) {
// oil-gas: we require SGAS
if (!deck.hasField("SGAS")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS keyword in 2-phase init");
}
const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
for (int c = 0; c < number_of_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
s[2*c + gpos] = sg_deck[c_deck];
s[2*c + opos] = 1.0 - sg_deck[c_deck];
p[c] = p_deck[c_deck];
}
} else {
// water-oil or water-gas: we require SWAT
if (!deck.hasField("SWAT")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SWAT keyword in 2-phase init");
}
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int nwpos = (wpos + 1) % 2;
for (int c = 0; c < number_of_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
s[2*c + wpos] = sw_deck[c_deck];
s[2*c + nwpos] = 1.0 - sw_deck[c_deck];
p[c] = p_deck[c_deck];
}
}
} else if (num_phases == 3) {
const bool has_swat_sgas = deck.hasField("SWAT") && deck.hasField("SGAS");
if (!has_swat_sgas) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init.");
}
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
for (int c = 0; c < number_of_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
s[3*c + wpos] = sw_deck[c_deck];
s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
s[3*c + gpos] = sg_deck[c_deck];
p[c] = p_deck[c_deck];
}
} else {
OPM_THROW(std::runtime_error, "initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases.");
}
} else {
OPM_THROW(std::runtime_error, "initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS.");
}
// Finally, init face pressures.
initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
begin_cell_centroids, state);
}
/// Initialize surface volume from pressure and saturation by z = As.
/// Here saturation is used as an intial guess for z in the
/// computation of A.
@ -975,60 +851,6 @@ namespace Opm
}
}
/// Initialize a blackoil state from input deck.
template <class Props, class State>
void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
const EclipseGridParser& deck,
const double gravity,
State& state)
{
initBlackoilStateFromDeck(grid.number_of_cells, grid.global_cell,
grid.number_of_faces, UgGridHelpers::faceCells(grid),
grid.face_centroids, grid.cell_centroids,grid.dimensions,
props, deck, gravity, state);
}
template <class FaceCells, class FCI, class CCI, class Props, class State>
void initBlackoilStateFromDeck(int number_of_cells,
const int* global_cell,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const Props& props,
const EclipseGridParser& deck,
const double gravity,
State& state)
{
initStateFromDeck(number_of_cells, global_cell, number_of_faces,
face_cells, begin_face_centroids, begin_cell_centroids,
dimensions, props, deck, gravity, state);
if (deck.hasField("RS")) {
const std::vector<double>& rs_deck = deck.getFloatingPointValue("RS");
for (int c = 0; c < number_of_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
state.gasoilratio()[c] = rs_deck[c_deck];
}
initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
computeSaturation(props,state);
} else if (deck.hasField("RV")){
const std::vector<double>& rv_deck = deck.getFloatingPointValue("RV");
for (int c = 0; c < number_of_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
state.rv()[c] = rv_deck[c_deck];
}
initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
computeSaturation(props,state);
}
else {
OPM_THROW(std::runtime_error, "Temporarily, we require the RS or the RV field.");
}
}
/// Initialize a blackoil state from input deck.
template <class Props, class State>
void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
@ -1084,7 +906,6 @@ namespace Opm
}
}
} // namespace Opm

44
opm/core/wells/WellCollection.cpp
View File

@ -79,50 +79,6 @@ namespace Opm
child->setParent(parent);
}
void WellCollection::addChild(const std::string& child_name,
const std::string& parent_name,
const EclipseGridParser& deck)
{
WellsGroupInterface* parent = findNode(parent_name);
if (!parent) {
roots_.push_back(createWellsGroup(parent_name, deck));
parent = roots_[roots_.size() - 1].get();
}
std::shared_ptr<WellsGroupInterface> child;
for (size_t i = 0; i < roots_.size(); ++i) {
if (roots_[i]->name() == child_name) {
child = roots_[i];
// We've found a new parent to the previously thought root, need to remove it
for(size_t j = i; j < roots_.size() - 1; ++j) {
roots_[j] = roots_[j+1];
}
roots_.resize(roots_.size()-1);
break;
}
}
if (!child.get()) {
child = createWellsGroup(child_name, deck);
}
WellsGroup* parent_as_group = static_cast<WellsGroup*> (parent);
if (!parent_as_group) {
OPM_THROW(std::runtime_error, "Trying to add child to group named " << parent_name << ", but it's not a group.");
}
parent_as_group->addChild(child);
if(child->isLeafNode()) {
leaf_nodes_.push_back(static_cast<WellNode*>(child.get()));
}
child->setParent(parent);
}
const std::vector<WellNode*>& WellCollection::getLeafNodes() const {
return leaf_nodes_;
}

14
opm/core/wells/WellCollection.hpp
View File

@ -16,9 +16,6 @@
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_WELLCOLLECTION_HPP
#define OPM_WELLCOLLECTION_HPP
@ -27,7 +24,6 @@
#include <opm/core/wells/WellsGroup.hpp>
#include <opm/core/grid.h>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
@ -47,16 +43,6 @@ namespace Opm
void addGroup(GroupConstPtr groupChild, std::string parent_name,
size_t timeStep, const PhaseUsage& phaseUsage);
/// Adds and creates if necessary the child to the collection
/// and appends it to parent's children. Also adds and creates the parent
/// if necessary.
/// \param[in] child name of child node
/// \param[in] parent name of parent node
/// \param[in] deck deck from which we will extract group control data
void addChild(const std::string& child,
const std::string& parent,
const EclipseGridParser& deck);
/// Adds the child to the collection
/// and appends it to parent's children.
/// \param[in] child the child node

95
opm/core/wells/WellsGroup.cpp
View File

@ -1045,101 +1045,6 @@ namespace Opm
}
} // anonymous namespace
std::shared_ptr<WellsGroupInterface> createWellsGroup(const std::string& name,
const EclipseGridParser& deck)
{
PhaseUsage phase_usage = phaseUsageFromDeck(deck);
std::shared_ptr<WellsGroupInterface> return_value;
// First we need to determine whether it's a group or just a well:
bool isWell = false;
if (deck.hasField("WELSPECS")) {
WELSPECS wspecs = deck.getWELSPECS();
for (size_t i = 0; i < wspecs.welspecs.size(); i++) {
if (wspecs.welspecs[i].name_ == name) {
isWell = true;
break;
}
}
}
// For now, assume that if it isn't a well, it's a group
if (isWell) {
ProductionSpecification production_specification;
InjectionSpecification injection_specification;
if (deck.hasField("WCONINJE")) {
WCONINJE wconinje = deck.getWCONINJE();
for (size_t i = 0; i < wconinje.wconinje.size(); i++) {
if (wconinje.wconinje[i].well_ == name) {
WconinjeLine line = wconinje.wconinje[i];
injection_specification.BHP_limit_ = line.BHP_limit_;
injection_specification.injector_type_ = toInjectorType(line.injector_type_);
injection_specification.control_mode_ = toInjectionControlMode(line.control_mode_);
injection_specification.surface_flow_max_rate_ = line.surface_flow_max_rate_;
injection_specification.reservoir_flow_max_rate_ = line.reservoir_flow_max_rate_;
production_specification.guide_rate_ = 0.0; // We know we're not a producer
}
}
}
if (deck.hasField("WCONPROD")) {
WCONPROD wconprod = deck.getWCONPROD();
for (size_t i = 0; i < wconprod.wconprod.size(); i++) {
if (wconprod.wconprod[i].well_ == name) {
WconprodLine line = wconprod.wconprod[i];
production_specification.BHP_limit_ = line.BHP_limit_;
production_specification.reservoir_flow_max_rate_ = line.reservoir_flow_max_rate_;
production_specification.oil_max_rate_ = line.oil_max_rate_;
production_specification.control_mode_ = toProductionControlMode(line.control_mode_);
production_specification.water_max_rate_ = line.water_max_rate_;
injection_specification.guide_rate_ = 0.0; // we know we're not an injector
}
}
}
return_value.reset(new WellNode(name, production_specification, injection_specification, phase_usage));
} else {
InjectionSpecification injection_specification;
if (deck.hasField("GCONINJE")) {
GCONINJE gconinje = deck.getGCONINJE();
for (size_t i = 0; i < gconinje.gconinje.size(); i++) {
if (gconinje.gconinje[i].group_ == name) {
GconinjeLine line = gconinje.gconinje[i];
injection_specification.injector_type_ = toInjectorType(line.injector_type_);
injection_specification.control_mode_ = toInjectionControlMode(line.control_mode_);
injection_specification.surface_flow_max_rate_ = line.surface_flow_max_rate_;
injection_specification.reservoir_flow_max_rate_ = line.resv_flow_max_rate_;
injection_specification.reinjection_fraction_target_ = line.reinjection_fraction_target_;
injection_specification.voidage_replacment_fraction_ = line.voidage_replacement_fraction_;
}
}
}
ProductionSpecification production_specification;
if (deck.hasField("GCONPROD")) {
std::cout << "Searching in gconprod " << std::endl;
std::cout << "name= " << name << std::endl;
GCONPROD gconprod = deck.getGCONPROD();
for (size_t i = 0; i < gconprod.gconprod.size(); i++) {
if (gconprod.gconprod[i].group_ == name) {
GconprodLine line = gconprod.gconprod[i];
production_specification.oil_max_rate_ = line.oil_max_rate_;
std::cout << "control_mode = " << line.control_mode_ << std::endl;
production_specification.control_mode_ = toProductionControlMode(line.control_mode_);
production_specification.water_max_rate_ = line.water_max_rate_;
production_specification.gas_max_rate_ = line.gas_max_rate_;
production_specification.liquid_max_rate_ = line.liquid_max_rate_;
production_specification.procedure_ = toProductionProcedure(line.procedure_);
production_specification.reservoir_flow_max_rate_ = line.resv_max_rate_;
}
}
}
return_value.reset(new WellsGroup(name, production_specification, injection_specification, phase_usage));
}
return return_value;
}
std::shared_ptr<WellsGroupInterface> createGroupWellsGroup(GroupConstPtr group, size_t timeStep, const PhaseUsage& phase_usage )
{
InjectionSpecification injection_specification;

7
opm/core/wells/WellsGroup.hpp
View File

@ -22,7 +22,6 @@
#include <opm/core/wells/InjectionSpecification.hpp>
#include <opm/core/wells/ProductionSpecification.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/grid.h>
#include <opm/core/props/BlackoilPhases.hpp>
@ -403,12 +402,6 @@ namespace Opm
bool shut_well_;
};
/// Creates the WellsGroupInterface for the given name
/// \param[in] name the name of the wells group.
/// \param[in] deck the deck from which to fetch information.
std::shared_ptr<WellsGroupInterface> createWellsGroup(const std::string& name,
const EclipseGridParser& deck);
/// Creates the WellsGroupInterface for the given well
/// \param[in] well the Well to construct object for
/// \param[in] timeStep the time step in question

1
opm/core/wells/WellsManager.cpp
View File

@ -21,7 +21,6 @@
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>

2
opm/core/wells/WellsManager.hpp
View File

@ -33,8 +33,6 @@ struct UnstructuredGrid;
namespace Opm
{
class EclipseGridParser;
struct WellData
{
WellType type;

12
tests/test_equil.cpp
View File

@ -337,7 +337,7 @@ BOOST_AUTO_TEST_CASE (DeckAllDead)
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("deadfluids.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, *grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, *grid, 10.0);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, *grid, 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells);
@ -416,7 +416,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
Opm::DeckConstPtr deck = parser->parseFile("capillary.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 10.0);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -455,7 +455,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
Opm::DeckConstPtr deck = parser->parseFile("capillary_overlap.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -516,7 +516,7 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
Opm::DeckConstPtr deck = parser->parseFile("equil_liveoil.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -594,7 +594,7 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
Opm::DeckConstPtr deck = parser->parseFile("equil_livegas.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -675,7 +675,7 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
Opm::DeckConstPtr deck = parser->parseFile("equil_rsvd_and_rvvd.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);