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

Refactored parts needed for Blackoil in autodiff to get rid of UG dependency.

Browse Source 
This patch refactors (hopefully) all parts of opm-core that are needed
by the fully implicite black oil solver in opm-autodiff and that inherently
relied on UnstructuredGrid.

We added a new simple grid interface consisting out of free functions
that will allow us to use CpGrid without copying it to an UnstructuredGrid
by the means of the GridAdapter. Using this interface we have add methods that
allow specifying the grid information (global_cell, cartdims, etc.) wherever
possible to prevent introducing grid parameters for the type of the grid.
Unfortunately this was not possible everywhere.
This commit is contained in:
Markus Blatt 2014-02-17 13:23:01 +01:00
parent 3d215691da
commit d5f470cb68
15 changed files with 667 additions and 288 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -28,83 +28,18 @@ namespace Opm
const UnstructuredGrid& grid,
bool init_rock)
{
if (init_rock){
rock_.init(deck, grid);
}
pvt_.init(deck, 200);
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, 200);
if (pvt_.numPhases() != satprops_->numPhases()) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
}
init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids,
grid.dimensions, init_rock);
}
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,
bool init_rock)
{
if(init_rock){
rock_.init(deck, grid);
}
const int pvt_samples = param.getDefault("pvt_tab_size", 200);
pvt_.init(deck, pvt_samples);
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", 200);
std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
if (deck.hasField("ENDSCALE") && threephase_model != "simple") {
OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'simple' model only.");
}
if (sat_samples > 1) {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegUniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else {
OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
}
} else {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else {
OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
}
}
if (pvt_.numPhases() != satprops_->numPhases()) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
}
init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids,
grid.dimensions, param, init_rock);
}
BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck()

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

@ -184,6 +184,24 @@ namespace Opm
double* smax) const;
private:
template<class T>
void init(const EclipseGridParser& deck,
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(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);
RockFromDeck rock_;
BlackoilPvtProperties pvt_;
std::unique_ptr<SaturationPropsInterface> satprops_;
@ -197,5 +215,6 @@ namespace Opm
} // namespace Opm
#include "BlackoilPropertiesFromDeck_impl.hpp"
#endif // OPM_BLACKOILPROPERTIESFROMDECK_HEADER_INCLUDED

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

@ -58,22 +58,34 @@ namespace Opm
void RockFromDeck::init(const EclipseGridParser& deck,
const UnstructuredGrid& grid)
{
assignPorosity(deck, grid);
permfield_valid_.assign(grid.number_of_cells, false);
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, grid, perm_threshold);
assignPermeability(deck, number_of_cells, global_cell, cart_dims, perm_threshold);
}
void RockFromDeck::assignPorosity(const EclipseGridParser& parser,
const UnstructuredGrid& grid)
int number_of_cells, const int* global_cell)
{
porosity_.assign(grid.number_of_cells, 1.0);
const int* gc = grid.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 = (gc == NULL) ? c : gc[c];
const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
porosity_[c] = poro[deck_pos];
}
}
@ -82,16 +94,17 @@ namespace Opm
void RockFromDeck::assignPermeability(const EclipseGridParser& parser,
const UnstructuredGrid& grid,
int number_of_cells,
const int* global_cell,
const int* cartdims,
double perm_threshold)
{
const int dim = 3;
const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2];
const int nc = grid.number_of_cells;
const int num_global_cells = cartdims[0]*cartdims[1]*cartdims[2];
assert(num_global_cells > 0);
permeability_.assign(dim * dim * nc, 0.0);
permeability_.assign(dim * dim * number_of_cells, 0.0);
std::vector<const std::vector<double>*> tensor;
tensor.reserve(10);
@ -113,13 +126,12 @@ namespace Opm
// chosen) default value...
//
if (tensor.size() > 1) {
const int* gc = grid.global_cell;
int off = 0;
for (int c = 0; c < nc; ++c, off += dim*dim) {
for (int c = 0; c < number_of_cells; ++c, off += dim*dim) {
// SharedPermTensor K(dim, dim, &permeability_[off]);
int kix = 0;
const int glob = (gc == NULL) ? c : gc[c];
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) {

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

@ -43,6 +43,16 @@ namespace Opm
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);
/// \return D, the number of spatial dimensions. Always 3 for deck input.
int numDimensions() const
{
@ -71,9 +81,12 @@ namespace Opm
private:
void assignPorosity(const EclipseGridParser& parser,
const UnstructuredGrid& grid);
int number_of_cells,
const int* global_cell);
void assignPermeability(const EclipseGridParser& parser,
const UnstructuredGrid& grid,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
const double perm_threshold);
std::vector<double> porosity_;

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

@ -60,6 +60,28 @@ namespace Opm
const UnstructuredGrid& grid,
const int samples);
/// Initialize from deck and grid.
/// \param[in] deck 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);
/// \return P, the number of phases.
int numPhases() const;
@ -123,11 +145,14 @@ namespace Opm
typedef SatFuncSet Funcs;
const Funcs& funcForCell(const int cell) const;
template<class T>
void initEPS(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const std::string& keyword,
std::vector<double>& scaleparam);
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam);
void relpermEPS(const double *s, const int cell, double *kr, double *dkrds= 0) const;
};

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

@ -46,6 +46,19 @@ namespace Opm
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>
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);
@ -61,11 +74,9 @@ namespace Opm
if (deck.hasField("SATNUM")) {
const std::vector<int>& satnum = deck.getIntegerValue("SATNUM");
satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
const int num_cells = grid.number_of_cells;
cell_to_func_.resize(num_cells);
const int* gc = grid.global_cell;
for (int cell = 0; cell < num_cells; ++cell) {
const int deck_pos = (gc == NULL) ? cell : gc[cell];
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;
}
}
@ -118,14 +129,22 @@ namespace Opm
}
}
do_eps_ = true;
initEPS(deck, grid, std::string("SWCR"), eps_.swcr_);
initEPS(deck, grid, std::string("SWL"), eps_.swl_);
initEPS(deck, grid, std::string("SWU"), eps_.swu_);
initEPS(deck, grid, std::string("SOWCR"), eps_.sowcr_);
initEPS(deck, grid, std::string("KRW"), eps_.krw_);
initEPS(deck, grid, std::string("KRWR"), eps_.krwr_);
initEPS(deck, grid, std::string("KRO"), eps_.kro_);
initEPS(deck, grid, std::string("KRORW"), eps_.krorw_);
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions,
std::string("SWCR"), eps_.swcr_);
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions,
std::string("SWL"), eps_.swl_);
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions,
std::string("SWU"), eps_.swu_);
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions,
std::string("SOWCR"), eps_.sowcr_);
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions,
std::string("KRW"), eps_.krw_);
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions,
std::string("KRWR"), eps_.krwr_);
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions,
std::string("KRO"), eps_.kro_);
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions,
std::string("KRORW"), eps_.krorw_);
}
}
@ -253,10 +272,45 @@ namespace Opm
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
namespace
{
template<class T>
const T* increment(T* cc, int i, int dim)
{
return cc+(i*dim);
}
template<class T>
T increment(const T& t, int i, int)
{
return t+i;
}
template<class T>
double getCoordinate(T* cc, int i)
{
return cc[i];
}
template<class T>
double getCoordinate(T t, int i)
{
return t->center()[i];
}
}
// Initialize saturation scaling parameter
template <class SatFuncSet>
template <class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam)
{
@ -273,29 +327,29 @@ namespace Opm
if (keyword == std::string("SWL")) {
if (useKeyword || deck.getENPTVD().mask_[0]) {
itab = 1;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
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")) {
if (useKeyword || deck.getENPTVD().mask_[1]) {
itab = 2;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
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")) {
if (useKeyword || deck.getENPTVD().mask_[2]) {
itab = 3;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
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("SOWCR")) {
if (useKeyword || deck.getENPTVD().mask_[3]) {
itab = 4;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sowcr_;
}
}else {
@ -308,29 +362,29 @@ namespace Opm
if (keyword == std::string("KRW")) {
if (useKeyword || deck.getENKRVD().mask_[0]) {
itab = 1;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwmax_;
}
} else if (keyword == std::string("KRO")) {
if (useKeyword || deck.getENKRVD().mask_[1]) {
itab = 2;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
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")) {
if (useKeyword || deck.getENKRVD().mask_[2]) {
itab = 3;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwr_;
}
} else if (keyword == std::string("KRORW")) {
if (useKeyword || deck.getENKRVD().mask_[3]) {
itab = 4;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorw_;
}
} else {
@ -346,10 +400,9 @@ namespace Opm
} else if (useKeyword) {
// Keyword values from deck
std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
const int* gc = grid.global_cell;
const std::vector<double>& val = deck.getFloatingPointValue(keyword);
for (int c = 0; c < int(scaleparam.size()); ++c) {
const int deck_pos = (gc == NULL) ? c : gc[c];
const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
scaleparam[c] = val[deck_pos];
}
} else {
@ -358,14 +411,13 @@ namespace Opm
deck.getENPTVD().write(std::cout);
else
deck.getENKRVD().write(std::cout);
const double* cc = grid.cell_centroids;
const int dim = grid.dimensions;
for (int cell = 0; cell < grid.number_of_cells; ++cell) {
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 = cc[dim*cell+dim-1];
double zc = getCoordinate(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);
}

13
opm/core/simulator/BlackoilState.cpp
View File

@ -4,13 +4,18 @@
using namespace Opm;
void
BlackoilState::init(const UnstructuredGrid& g, int num_phases) {
SimulatorState::init(g, num_phases);
gor_.resize(g.number_of_cells, 0.) ;
rv_.resize(g.number_of_cells,0.);
BlackoilState::init(int number_of_cells, int number_of_phases, int num_phases) {
SimulatorState::init(number_of_cells, number_of_phases, num_phases);
gor_.resize(number_of_cells, 0.) ;
rv_.resize(number_of_cells,0.);
// surfvol_ intentionally empty, left to initBlackoilSurfvol
}
void
BlackoilState::init(const UnstructuredGrid& g, int num_phases)
{
init(g.number_of_cells, g.number_of_faces, num_phases);
}
/// Set the first saturation to either its min or max value in
/// the indicated cells. The second saturation value s2 is set
/// to (1.0 - s1) for each cell. Any further saturation values

4
opm/core/simulator/BlackoilState.hpp
View File

@ -32,7 +32,9 @@ namespace Opm
class BlackoilState : public SimulatorState
{
public:
virtual void init(const UnstructuredGrid& g, int num_phases);
virtual void init(const UnstructuredGrid& grid, int num_phases);
virtual void init(int number_of_cells, int number_of_faces, int num_phases);
/// Set the first saturation to either its min or max value in
/// the indicated cells. The second saturation value s2 is set

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

@ -65,6 +65,44 @@ namespace Opm
const double gravity,
State& state);
/// Initialize a two-phase state from parameters.
/// The following parameters are accepted (defaults):
/// - convection_testcase (false) -- Water in the 'left' part of the grid.
/// - ref_pressure (100) -- Initial pressure in bar for all cells
/// (if convection_testcase is true),
/// or pressure at woc depth.
/// - segregation_testcase (false) -- Water above the woc instead of below.
/// - water_oil_contact (none) -- Depth of water-oil contact (woc).
/// - init_saturation (none) -- Initial water saturation for all cells.
///
/// If convection_testcase is true, the saturation is initialised
/// as indicated, and pressure is initialised to a constant value
/// ('ref_pressure').
/// If segregation_testcase is true, the saturation is initialised
/// as indicated, and pressure is initialised hydrostatically.
/// Otherwise we have 3 cases:
/// 1. If 'water_oil_contact' is given, saturation is initialised
/// accordingly.
/// 2. If 'water_oil_contact' is not given, but 'init_saturation'
/// is given, water saturation is set to that value everywhere.
/// 3. If neither are given, water saturation is set to minimum.
///
/// In all three cases, pressure is initialised hydrostatically.
/// In case 2) and 3), the depth of the first cell is used as reference depth.
template <class FaceCells, class CCI, class FCI, class State>
void initStateBasic(int number_of_cells,
const int* global_cell,
const int* cartdims,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const IncompPropertiesInterface& props,
const parameter::ParameterGroup& param,
const double gravity,
State& state);
/// Initialize a blackoil state from parameters.
/// The following parameters are accepted (defaults):
/// - convection_testcase (false) -- Water in the 'left' part of the grid.
@ -88,6 +126,36 @@ namespace Opm
const double gravity,
State& state);
/// Initialize a blackoil state from parameters.
/// The following parameters are accepted (defaults):
/// - convection_testcase (false) -- Water in the 'left' part of the grid.
/// - ref_pressure (100) -- Initial pressure in bar for all cells
/// (if convection_testcase is true),
/// or pressure at woc depth.
/// - water_oil_contact (none) -- Depth of water-oil contact (woc).
/// If convection_testcase is true, the saturation is initialised
/// as indicated, and pressure is initialised to a constant value
/// ('ref_pressure').
/// Otherwise we have 2 cases:
/// 1. If 'water_oil_contact' is given, saturation is initialised
/// accordingly.
/// 2. Water saturation is set to minimum.
/// In both cases, pressure is initialised hydrostatically.
/// In case 2., the depth of the first cell is used as reference depth.
template <class FaceCells, class FCI, class CCI, class State>
void initStateBasic(int number_of_cells,
const int* global_cell,
const int* cartdims,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const BlackoilPropertiesInterface& props,
const parameter::ParameterGroup& param,
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,
@ -102,6 +170,27 @@ namespace Opm
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,
const int* cartdims,
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 two-phase water-oil blackoil state from input deck.
/// If EQUIL is present:
/// - saturation is set according to the water-oil contact,
@ -117,7 +206,27 @@ namespace Opm
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,
const int* cartdims,
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);
} // namespace Opm
#include <opm/core/simulator/initState_impl.hpp>

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

@ -23,6 +23,7 @@
#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>
@ -47,17 +48,21 @@ namespace Opm
// Find the cells that are below and above a depth.
// TODO: add 'anitialiasing', obtaining a more precise split
// by f. ex. subdividing cells cut by the split depths.
void cellsBelowAbove(const UnstructuredGrid& grid,
template<class T>
void cellsBelowAbove(int number_of_cells,
T begin_cell_centroids,
int dimension,
const double depth,
std::vector<int>& below,
std::vector<int>& above)
{
const int num_cells = grid.number_of_cells;
below.reserve(num_cells);
above.reserve(num_cells);
const int dim = grid.dimensions;
for (int c = 0; c < num_cells; ++c) {
const double z = grid.cell_centroids[dim*c + dim - 1];
below.reserve(number_of_cells);
above.reserve(number_of_cells);
for (int c = 0; c < number_of_cells; ++c) {
const double z = UgGridHelpers
::getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c,
dimension),
dimension-1);
if (z > depth) {
below.push_back(c);
} else {
@ -74,8 +79,10 @@ namespace Opm
// Initialize saturations so that there is water below woc,
// and oil above.
// If invert is true, water is instead above, oil below.
template <class Props, class State>
void initWaterOilContact(const UnstructuredGrid& grid,
template <class T, class Props, class State>
void initWaterOilContact(int number_of_cells,
T begin_cell_centroids,
int dimensions,
const Props& props,
const double woc,
const WaterInit waterinit,
@ -91,10 +98,12 @@ namespace Opm
// }
switch (waterinit) {
case WaterBelow:
cellsBelowAbove(grid, woc, water, oil);
cellsBelowAbove(number_of_cells, begin_cell_centroids, dimensions,
woc, water, oil);
break;
case WaterAbove:
cellsBelowAbove(grid, woc, oil, water);
cellsBelowAbove(number_of_cells, begin_cell_centroids, dimensions,
woc, oil, water);
}
// Set saturations.
state.setFirstSat(oil, props, State::MinSat);
@ -113,8 +122,10 @@ namespace Opm
// Note that by manipulating the densities parameter,
// it is possible to initialise with water on top instead of
// on the bottom etc.
template <class State>
void initHydrostaticPressure(const UnstructuredGrid& grid,
template <class T, class State>
void initHydrostaticPressure(int number_of_cells,
T begin_cell_centroids,
int dimensions,
const double* densities,
const double woc,
const double gravity,
@ -123,14 +134,14 @@ namespace Opm
State& state)
{
std::vector<double>& p = state.pressure();
const int num_cells = grid.number_of_cells;
const int dim = grid.dimensions;
// Compute pressure at woc
const double rho_datum = datum_z > woc ? densities[0] : densities[1];
const double woc_p = datum_p + (woc - datum_z)*gravity*rho_datum;
for (int c = 0; c < num_cells; ++c) {
for (int c = 0; c < number_of_cells; ++c) {
// Compute pressure as delta from woc pressure.
const double z = grid.cell_centroids[dim*c + dim - 1];
const double z = UgGridHelpers::
getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions),
dimensions-1);
const double rho = z > woc ? densities[0] : densities[1];
p[c] = woc_p + (z - woc)*gravity*rho;
}
@ -139,8 +150,10 @@ namespace Opm
// Facade to initHydrostaticPressure() taking a property object,
// for similarity to initHydrostaticPressure() for compressible fluids.
template <class State>
void initHydrostaticPressure(const UnstructuredGrid& grid,
template <class T, class State>
void initHydrostaticPressure(int number_of_cells,
T begin_cell_centroids,
int dimensions,
const IncompPropertiesInterface& props,
const double woc,
const double gravity,
@ -149,7 +162,8 @@ namespace Opm
State& state)
{
const double* densities = props.density();
initHydrostaticPressure(grid, densities, woc, gravity, datum_z, datum_p, state);
initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
densities, woc, gravity, datum_z, datum_p, state);
}
@ -181,8 +195,10 @@ namespace Opm
// where rho is the oil density above the woc, water density
// below woc. Note that even if there is (immobile) water in
// the oil zone, it does not contribute to the pressure there.
template <class State>
void initHydrostaticPressure(const UnstructuredGrid& grid,
template <class T, class State>
void initHydrostaticPressure(int number_of_cells,
T begin_cell_centroids,
int dimensions,
const BlackoilPropertiesInterface& props,
const double woc,
const double gravity,
@ -193,11 +209,11 @@ namespace Opm
assert(props.numPhases() == 2);
// Obtain max and min z for which we will need to compute p.
const int num_cells = grid.number_of_cells;
const int dim = grid.dimensions;
double zlim[2] = { 1e100, -1e100 };
for (int c = 0; c < num_cells; ++c) {
const double z = grid.cell_centroids[dim*c + dim - 1];
for (int c = 0; c < number_of_cells; ++c) {
const double z = UgGridHelpers::
getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions),
dimensions-1);;
zlim[0] = std::min(zlim[0], z);
zlim[1] = std::max(zlim[1], z);
}
@ -266,30 +282,42 @@ namespace Opm
// Evaluate pressure at each cell centroid.
std::vector<double>& p = state.pressure();
for (int c = 0; c < num_cells; ++c) {
const double z = grid.cell_centroids[dim*c + dim - 1];
for (int c = 0; c < number_of_cells; ++c) {
const double z = UgGridHelpers::
getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions),
dimensions-1);
p[c] = press(z);
}
}
// Initialize face pressures to distance-weighted average of adjacent cell pressures.
template <class State>
void initFacePressure(const UnstructuredGrid& grid,
template <class State, class FaceCells, class FCI, class CCI>
void initFacePressure(int dimensions,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
State& state)
{
const int dim = grid.dimensions;
const std::vector<double>& cp = state.pressure();
std::vector<double>& fp = state.facepressure();
for (int f = 0; f < grid.number_of_faces; ++f) {
for (int f = 0; f < number_of_faces; ++f) {
double dist[2] = { 0.0, 0.0 };
double press[2] = { 0.0, 0.0 };
int bdy_idx = -1;
for (int j = 0; j < 2; ++j) {
const int c = grid.face_cells[2*f + j];
const int c = face_cells(f, j);
if (c >= 0) {
dist[j] = 0.0;
for (int dd = 0; dd < dim; ++dd) {
double diff = grid.face_centroids[dim*f + dd] - grid.cell_centroids[dim*c + dd];
for (int dd = 0; dd < dimensions; ++dd) {
double diff = UgGridHelpers::
getCoordinate(UgGridHelpers::increment(begin_face_centroids, f,
dimensions),
dd)
- UgGridHelpers::
getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c,
dimensions),
dd);
dist[j] += diff*diff;
}
dist[j] = std::sqrt(dist[j]);
@ -318,11 +346,32 @@ namespace Opm
const double gravity,
State& state)
{
initStateBasic(grid.number_of_cells, grid.global_cell, grid.cartdims,
grid.number_of_faces, UgGridHelpers::faceCells(grid),
grid.face_centroids, grid.cell_centroids,
grid.dimensions, props, param,
gravity, state);
}
template <class FaceCells, class CCI, class FCI, class State>
void initStateBasic(int number_of_cells,
const int* global_cell,
const int* cartdims,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const IncompPropertiesInterface& props,
const parameter::ParameterGroup& param,
const double gravity,
State& state)
{
const int num_phases = props.numPhases();
if (num_phases != 2) {
OPM_THROW(std::runtime_error, "initStateTwophaseBasic(): currently handling only two-phase scenarios.");
}
state.init(grid, num_phases);
state.init(number_of_cells, number_of_faces, num_phases);
const int num_cells = props.numCells();
// By default: initialise water saturation to minimum everywhere.
std::vector<int> all_cells(num_cells);
@ -336,11 +385,9 @@ namespace Opm
// Initialise water saturation to max in the 'left' part.
std::vector<int> left_cells;
left_cells.reserve(num_cells/2);
const int *glob_cell = grid.global_cell;
const int* cd = grid.cartdims;
for (int cell = 0; cell < num_cells; ++cell) {
const int gc = glob_cell == 0 ? cell : glob_cell[cell];
bool left = (gc % cd[0]) < cd[0]/2;
const int gc = global_cell == 0 ? cell : global_cell[cell];
bool left = (gc % cartdims[0]) < cartdims[0]/2;
if (left) {
left_cells.push_back(cell);
}
@ -353,30 +400,34 @@ namespace Opm
if (gravity == 0.0) {
std::cout << "**** Warning: running gravity segregation scenario, but gravity is zero." << std::endl;
}
if (grid.cartdims[2] <= 1) {
if (cartdims[2] <= 1) {
std::cout << "**** Warning: running gravity segregation scenario, which expects nz > 1." << std::endl;
}
// Initialise water saturation to max *above* water-oil contact.
const double woc = param.get<double>("water_oil_contact");
initWaterOilContact(grid, props, woc, WaterAbove, state);
initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
props, woc, WaterAbove, state);
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
double dens[2] = { props.density()[1], props.density()[0] };
initHydrostaticPressure(grid, dens, woc, gravity, woc, ref_p, state);
initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
dens, woc, gravity, woc, ref_p, state);
} else if (param.has("water_oil_contact")) {
// Warn against error-prone usage.
if (gravity == 0.0) {
std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
}
if (grid.cartdims[2] <= 1) {
if (cartdims[2] <= 1) {
std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl;
}
// Initialise water saturation to max below water-oil contact.
const double woc = param.get<double>("water_oil_contact");
initWaterOilContact(grid, props, woc, WaterBelow, state);
initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
props, woc, WaterBelow, state);
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
initHydrostaticPressure(grid, props.density(), woc, gravity, woc, ref_p, state);
initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
props.density(), woc, gravity, woc, ref_p, state);
} else if (param.has("init_saturation")) {
// Initialise water saturation to init_saturation parameter.
const double init_saturation = param.get<double>("init_saturation");
@ -388,20 +439,23 @@ namespace Opm
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
const double rho = props.density()[0]*init_saturation + props.density()[1]*(1.0 - init_saturation);
const double dens[2] = { rho, rho };
const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1];
initHydrostaticPressure(grid, dens, ref_z, gravity, ref_z, ref_p, state);
const double ref_z = getCoordinate(begin_cell_centroids, dimensions - 1);
initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
dens, ref_z, gravity, ref_z, ref_p, state);
} else {
// Use default: water saturation is minimum everywhere.
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
const double rho = props.density()[1];
const double dens[2] = { rho, rho };
const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1];
initHydrostaticPressure(grid, dens, ref_z, gravity, ref_z, ref_p, state);
const double ref_z = getCoordinate(begin_cell_centroids, dimensions - 1);
initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
dens, ref_z, gravity, ref_z, ref_p, state);
}
// Finally, init face pressures.
initFacePressure(grid, state);
initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
begin_cell_centroids, state);
}
@ -412,13 +466,33 @@ namespace Opm
const parameter::ParameterGroup& param,
const double gravity,
State& state)
{
initStateBasic(grid.number_of_cells, grid.global_cell, grid.cartdims,
grid.number_of_faces, UgGridHelpers::faceCells(grid),
grid.face_centroids, grid.cell_centroids, grid.dimensions,
props, param, gravity, state);
}
template <class FaceCells, class FCI, class CCI, class State>
void initStateBasic(int number_of_cells,
const int* global_cell,
const int* cartdims,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const BlackoilPropertiesInterface& props,
const parameter::ParameterGroup& param,
const double gravity,
State& state)
{
// TODO: Refactor to exploit similarity with IncompProp* case.
const int num_phases = props.numPhases();
if (num_phases != 2) {
OPM_THROW(std::runtime_error, "initStateTwophaseBasic(): currently handling only two-phase scenarios.");
}
state.init(grid, num_phases);
state.init(number_of_cells, number_of_faces, num_phases);
const int num_cells = props.numCells();
// By default: initialise water saturation to minimum everywhere.
std::vector<int> all_cells(num_cells);
@ -431,11 +505,9 @@ namespace Opm
// Initialise water saturation to max in the 'left' part.
std::vector<int> left_cells;
left_cells.reserve(num_cells/2);
const int *glob_cell = grid.global_cell;
const int* cd = grid.cartdims;
for (int cell = 0; cell < num_cells; ++cell) {
const int gc = glob_cell == 0 ? cell : glob_cell[cell];
bool left = (gc % cd[0]) < cd[0]/2;
const int gc = global_cell == 0 ? cell : global_cell[cell];
bool left = (gc % cartdims[0]) < cartdims[0]/2;
if (left) {
left_cells.push_back(cell);
}
@ -448,26 +520,30 @@ namespace Opm
if (gravity == 0.0) {
std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
}
if (grid.cartdims[2] <= 1) {
if (cartdims[2] <= 1) {
std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl;
}
// Initialise water saturation to max below water-oil contact.
const double woc = param.get<double>("water_oil_contact");
initWaterOilContact(grid, props, woc, WaterBelow, state);
initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
props, woc, WaterBelow, state);
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
initHydrostaticPressure(grid, props, woc, gravity, woc, ref_p, state);
initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
props, woc, gravity, woc, ref_p, state);
} else {
// Use default: water saturation is minimum everywhere.
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1];
const double ref_z = UgGridHelpers::getCoordinate(begin_cell_centroids, dimensions - 1);
const double woc = -1e100;
initHydrostaticPressure(grid, props, woc, gravity, ref_z, ref_p, state);
initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
props, woc, gravity, ref_z, ref_p, state);
}
// Finally, init face pressures.
initFacePressure(grid, state);
initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
begin_cell_centroids, state);
}
@ -478,6 +554,26 @@ namespace Opm
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);
@ -485,7 +581,7 @@ namespace Opm
OPM_THROW(std::runtime_error, "initStateFromDeck(): user specified property object with " << num_phases << " phases, "
"found " << pu.num_phases << " phases in deck.");
}
state.init(grid, num_phases);
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.");
@ -499,17 +595,18 @@ namespace Opm
OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet.");
}
const double woc = equil.equil[0].water_oil_contact_depth_;
initWaterOilContact(grid, props, woc, WaterBelow, state);
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(grid, props, woc, gravity, datum_z, datum_p, state);
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");
const int num_cells = grid.number_of_cells;
if (num_phases == 2) {
if (!pu.phase_used[BlackoilPhases::Aqua]) {
// oil-gas: we require SGAS
@ -519,8 +616,8 @@ namespace Opm
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 < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
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];
@ -533,8 +630,8 @@ namespace Opm
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 < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
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];
@ -550,8 +647,8 @@ namespace Opm
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 < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
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];
@ -565,7 +662,8 @@ namespace Opm
}
// Finally, init face pressures.
initFacePressure(grid, state);
initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
begin_cell_centroids, state);
}
@ -579,10 +677,18 @@ namespace Opm
void initBlackoilSurfvol(const UnstructuredGrid& grid,
const Props& props,
State& state)
{
initBlackoilSurfvol(grid.number_of_cells, props, state);
}
template <class Props, class State>
void initBlackoilSurfvol(int number_of_cells,
const Props& props,
State& state)
{
state.surfacevol() = state.saturation();
const int np = props.numPhases();
const int nc = grid.number_of_cells;
const int nc = number_of_cells;
std::vector<double> allA(nc*np*np);
std::vector<int> allcells(nc);
for (int c = 0; c < nc; ++c) {
@ -616,6 +722,15 @@ namespace Opm
const Props& props,
State& state)
{
initBlackoilSurfvolUsingRSorRV(grid.number_of_cells, props, state);
}
template <class Props, class State>
void initBlackoilSurfvolUsingRSorRV(int number_of_cells,
const Props& props,
State& state)
{
if (props.numPhases() != 3) {
OPM_THROW(std::runtime_error, "initBlackoilSurfvol() is only supported in three-phase simulations.");
}
@ -627,27 +742,26 @@ namespace Opm
const PhaseUsage pu = props.phaseUsage();
const int np = props.numPhases();
const int nc = grid.number_of_cells;
std::vector<double> allA_a(nc*np*np);
std::vector<double> allA_l(nc*np*np);
std::vector<double> allA_v(nc*np*np);
std::vector<double> allA_a(number_of_cells*np*np);
std::vector<double> allA_l(number_of_cells*np*np);
std::vector<double> allA_v(number_of_cells*np*np);
std::vector<int> allcells(nc);
std::vector<double> z_init(nc*np);
std::vector<int> allcells(number_of_cells);
std::vector<double> z_init(number_of_cells*np);
std::fill(z_init.begin(),z_init.end(),0.0);
for (int c = 0; c < nc; ++c) {
for (int c = 0; c < number_of_cells; ++c) {
allcells[c] = c;
}
std::vector<double> capPressures(nc*np);
props.capPress(nc,&state.saturation()[0],&allcells[0],&capPressures[0],NULL);
std::vector<double> capPressures(number_of_cells*np);
props.capPress(number_of_cells,&state.saturation()[0],&allcells[0],&capPressures[0],NULL);
std::vector<double> Pw(nc);
std::vector<double> Pg(nc);
std::vector<double> Pw(number_of_cells);
std::vector<double> Pg(number_of_cells);
for (int c = 0; c < nc; ++c){
for (int c = 0; c < number_of_cells; ++c){
Pw[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Aqua];
Pg[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Vapour];
}
@ -657,7 +771,7 @@ namespace Opm
// Water phase
if(pu.phase_used[BlackoilPhases::Aqua])
for (int c = 0; c < nc ; ++c){
for (int c = 0; c < number_of_cells ; ++c){
for (int p = 0; p < np ; ++p){
if (p == BlackoilPhases::Aqua)
z_tmp = 1;
@ -667,11 +781,11 @@ namespace Opm
z_init[c*np + p] = z_tmp;
}
}
props.matrix(nc, &state.pressure()[0], &z_init[0], &allcells[0], &allA_a[0], 0);
props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_a[0], 0);
// Liquid phase
if(pu.phase_used[BlackoilPhases::Liquid]){
for (int c = 0; c < nc ; ++c){
for (int c = 0; c < number_of_cells ; ++c){
for (int p = 0; p < np ; ++p){
if(p == BlackoilPhases::Vapour){
if(state.saturation()[np*c + p] > 0)
@ -689,10 +803,10 @@ namespace Opm
}
}
}
props.matrix(nc, &state.pressure()[0], &z_init[0], &allcells[0], &allA_l[0], 0);
props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_l[0], 0);
if(pu.phase_used[BlackoilPhases::Vapour]){
for (int c = 0; c < nc ; ++c){
for (int c = 0; c < number_of_cells ; ++c){
for (int p = 0; p < np ; ++p){
if(p == BlackoilPhases::Liquid){
if(state.saturation()[np*c + p] > 0)
@ -710,9 +824,9 @@ namespace Opm
}
}
}
props.matrix(nc, &state.pressure()[0], &z_init[0], &allcells[0], &allA_v[0], 0);
props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_v[0], 0);
for (int c = 0; c < nc; ++c) {
for (int c = 0; c < number_of_cells; ++c) {
// Using z = As
double* z = &state.surfacevol()[c*np];
const double* s = &state.saturation()[c*np];
@ -744,24 +858,44 @@ namespace Opm
const double gravity,
State& state)
{
initStateFromDeck(grid, props, deck, gravity, state);
initBlackoilStateFromDeck(grid.number_of_cells, grid.global_cell, grid.cartdims,
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,
const int* cartdims,
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");
const int num_cells = grid.number_of_cells;
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
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(grid, props, state);
initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
computeSaturation(props,state);
} else if (deck.hasField("RV")){
const std::vector<double>& rv_deck = deck.getFloatingPointValue("RV");
const int num_cells = grid.number_of_cells;
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
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(grid, props, state);
initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
computeSaturation(props,state);
}

28
opm/core/utility/miscUtilities.cpp
View File

@ -401,24 +401,16 @@ namespace Opm
const std::vector<double>& face_flux,
std::vector<double>& cell_velocity)
{
const int dim = grid.dimensions;
cell_velocity.clear();
cell_velocity.resize(grid.number_of_cells*dim, 0.0);
for (int face = 0; face < grid.number_of_faces; ++face) {
int c[2] = { grid.face_cells[2*face], grid.face_cells[2*face + 1] };
const double* fc = &grid.face_centroids[face*dim];
double flux = face_flux[face];
for (int i = 0; i < 2; ++i) {
if (c[i] >= 0) {
const double* cc = &grid.cell_centroids[c[i]*dim];
for (int d = 0; d < dim; ++d) {
double v_contrib = fc[d] - cc[d];
v_contrib *= flux/grid.cell_volumes[c[i]];
cell_velocity[c[i]*dim + d] += (i == 0) ? v_contrib : -v_contrib;
}
}
}
}
estimateCellVelocity(grid.number_of_cells,
grid.number_of_faces,
grid.face_centroids,
UgGridHelpers::faceCells(grid),
grid.cell_centroids,
grid.cell_volumes,
grid.dimensions,
face_flux,
cell_velocity);
}
/// Extract a vector of water saturations from a vector of

21
opm/core/utility/miscUtilities.hpp
View File

@ -188,6 +188,26 @@ namespace Opm
const std::vector<double>& face_flux,
std::vector<double>& cell_velocity);
/// @brief Estimates a scalar cell velocity from face fluxes.
/// @param[in] number_of_cells The number of cells of the grid
/// @param[in] number_of_faces The number of cells of the grid
/// @param[in] begin_face_centroids Iterator pointing to first face centroid.
/// @param[in] face_cells Mapping from faces to connected cells.
/// @param[in] dimensions The dimensions of the grid.
/// @param[in] begin_cell_centroids Iterator pointing to first cell centroid.
/// @param[in] face_flux signed per-face fluxes
/// @param[out] cell_velocity the estimated velocities.
template<class CC, class FC, class FC1, class CV>
void estimateCellVelocity(int number_of_cells,
int number_of_faces,
FC begin_face_centroids,
FC1 face_cells,
CC begin_cell_centroids,
CV begin_cell_volumes,
int dimension,
const std::vector<double>& face_flux,
std::vector<double>& cell_velocity);
/// Extract a vector of water saturations from a vector of
/// interleaved water and oil saturations.
void toWaterSat(const std::vector<double>& sboth,
@ -320,4 +340,5 @@ namespace Opm
} // namespace Opm
#include "miscUtilities_impl.hpp"
#endif // OPM_MISCUTILITIES_HEADER_INCLUDED

41
opm/core/utility/miscUtilities_impl.hpp Normal file
View File

@ -0,0 +1,41 @@
#include <opm/core/grid/GridHelpers.hpp>
namespace Opm
{
/// @brief Estimates a scalar cell velocity from face fluxes.
/// @param[in] number_of_cells The number of cells of the grid
/// @param[in] begin_face_centroids Iterator pointing to first face centroid.
/// @param[in] face_cells Mapping from faces to connected cells.
/// @param[in] dimensions The dimensions of the grid.
/// @param[in] begin_cell_centroids Iterator pointing to first cell centroid.
/// @param[in] face_flux signed per-face fluxes
/// @param[out] cell_velocity the estimated velocities.
template<class CC, class FC, class FC1, class CV>
void estimateCellVelocity(int number_of_cells,
int number_of_faces,
FC begin_face_centroids,
FC1 face_cells,
CC begin_cell_centroids,
CV begin_cell_volumes,
int dimension,
const std::vector<double>& face_flux,
std::vector<double>& cell_velocity)
{
cell_velocity.clear();
cell_velocity.resize(number_of_cells*dimension, 0.0);
for (int face = 0; face < number_of_faces; ++face) {
int c[2] = { face_cells(face, 0), face_cells(face, 1) };
FC fc = UgGridHelpers::increment(begin_face_centroids, face, dimension);
double flux = face_flux[face];
for (int i = 0; i < 2; ++i) {
if (c[i] >= 0) {
CC cc = UgGridHelpers::increment(begin_cell_centroids, c[i], dimension);
for (int d = 0; d < dimension; ++d) {
double v_contrib = UgGridHelpers::getCoordinate(fc, d) - UgGridHelpers::getCoordinate(cc, d);
v_contrib *= flux/begin_cell_volumes[c[i]];
cell_velocity[c[i]*dimension + d] += (i == 0) ? v_contrib : -v_contrib;
}
}
}
}
}
}

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

@ -23,6 +23,7 @@
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/utility/ErrorMacros.hpp>
@ -182,13 +183,17 @@ namespace
{
using namespace std;
std::array<double, 3> cube;
int num_local_faces = grid.cell_facepos[cell + 1] - grid.cell_facepos[cell];
typedef Opm::UgGridHelpers::Cell2FacesTraits<UnstructuredGrid>::Type Cell2Faces;
typedef Cell2Faces::row_type FaceRow;
Cell2Faces c2f=Opm::UgGridHelpers::cell2Faces(grid);
FaceRow faces=c2f[cell];
int num_local_faces = faces.size();
vector<double> x(num_local_faces);
vector<double> y(num_local_faces);
vector<double> z(num_local_faces);
for (int lf=0; lf<num_local_faces; ++ lf) {
int face = grid.cell_faces[grid.cell_facepos[cell] + lf];
const double* centroid = &grid.face_centroids[grid.dimensions*face];
int face = faces[lf];
const double* centroid = Opm::UgGridHelpers::faceCentroid(grid, face);
x[lf] = centroid[0];
y[lf] = centroid[1];
z[lf] = centroid[2];
@ -263,8 +268,8 @@ namespace Opm
/// Construct from existing wells object.
WellsManager::WellsManager(struct Wells* W)
: w_(clone_wells(W))
WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
: w_(clone_wells(W)), checkCellExistence_(checkCellExistence)
{
}
@ -274,10 +279,11 @@ namespace Opm
const size_t timeStep,
const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,
const double* permeability)
: w_(0)
const double* permeability,
bool checkCellExistence)
: w_(0), checkCellExistence_(checkCellExistence)
{
if (grid.dimensions != 3) {
if (UgGridHelpers::dimensions(grid) != 3) {
OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional.");
}
@ -287,7 +293,8 @@ namespace Opm
}
std::map<int,int> cartesian_to_compressed;
setupCompressedToCartesian(grid, cartesian_to_compressed);
setupCompressedToCartesian(UgGridHelpers::globalCell(grid),
UgGridHelpers::numCells(grid), cartesian_to_compressed);
// Obtain phase usage data.
PhaseUsage pu = phaseUsageFromDeck(eclipseState);
@ -417,10 +424,11 @@ namespace Opm
/// Construct wells from deck.
WellsManager::WellsManager(const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,
const double* permeability)
: w_(0)
const double* permeability,
bool checkCellExistence)
: w_(0), checkCellExistence_(checkCellExistence)
{
if (grid.dimensions != 3) {
if (UgGridHelpers::dimensions(grid) != 3) {
OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional.");
}
// NOTE: Implementation copied and modified from dune-porsol's class BlackoilWells.
@ -489,17 +497,17 @@ namespace Opm
// global_cell is a map from compressed cells to Cartesian grid cells.
// We must make the inverse lookup.
const int* global_cell = grid.global_cell;
const int* cpgdim = grid.cartdims;
const int* global_cell = UgGridHelpers::globalCell(grid);
const int* cpgdim = UgGridHelpers::cartDims(grid);
std::map<int,int> cartesian_to_compressed;
if (global_cell) {
for (int i = 0; i < grid.number_of_cells; ++i) {
for (int i = 0; i < UgGridHelpers::numCells(grid); ++i) {
cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
}
}
else {
for (int i = 0; i < grid.number_of_cells; ++i) {
for (int i = 0; i < UgGridHelpers::numCells(grid); ++i) {
cartesian_to_compressed.insert(std::make_pair(i, i));
}
}
@ -553,8 +561,10 @@ namespace Opm
std::map<int, int>::const_iterator cgit =
cartesian_to_compressed.find(cart_grid_indx);
if (cgit == cartesian_to_compressed.end()) {
OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << ix << ' ' << jy << ' '
<< kz << " not found in grid (well = " << name << ')');
if(checkCellExistence)
OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << ix << ' ' << jy << ' '
<< kz << " not found in grid (well = " << name << ')');
continue;
}
int cell = cgit->second;
PerfData pd;
@ -568,7 +578,8 @@ namespace Opm
OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
}
std::array<double, 3> cubical = getCubeDim(grid, cell);
const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell];
int dimensions = UgGridHelpers::dimensions(grid);
const double* cell_perm = &permeability[dimensions*dimensions*cell];
pd.well_index = computeWellIndex(radius, cubical, cell_perm,
compdat.compdat[kw].skin_factor_);
}
@ -586,7 +597,7 @@ namespace Opm
// Set up reference depths that were defaulted. Count perfs.
int num_perfs = 0;
assert(grid.dimensions == 3);
assert(UgGridHelpers::dimensions(grid) == 3);
for (int w = 0; w < num_wells; ++w) {
num_perfs += wellperf_data[w].size();
if (well_data[w].reference_bhp_depth < 0.0) {
@ -594,7 +605,8 @@ namespace Opm
double min_depth = 1e100;
int num_wperfs = wellperf_data[w].size();
for (int perf = 0; perf < num_wperfs; ++perf) {
double depth = grid.cell_centroids[3*wellperf_data[w][perf].cell + 2];
double depth = UgGridHelpers::
cellCentroidCoordinate(grid, wellperf_data[w][perf].cell, 2);
min_depth = std::min(min_depth, depth);
}
well_data[w].reference_bhp_depth = min_depth;
@ -1058,18 +1070,18 @@ namespace Opm
well_collection_.applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase);
}
void WellsManager::setupCompressedToCartesian(const UnstructuredGrid& grid, std::map<int,int>& cartesian_to_compressed ) {
void WellsManager::setupCompressedToCartesian(const int* global_cell, int number_of_cells,
std::map<int,int>& cartesian_to_compressed ) {
// global_cell is a map from compressed cells to Cartesian grid cells.
// We must make the inverse lookup.
const int* global_cell = grid.global_cell;
if (global_cell) {
for (int i = 0; i < grid.number_of_cells; ++i) {
for (int i = 0; i < number_of_cells; ++i) {
cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
}
}
else {
for (int i = 0; i < grid.number_of_cells; ++i) {
for (int i = 0; i < number_of_cells; ++i) {
cartesian_to_compressed.insert(std::make_pair(i, i));
}
}
@ -1118,12 +1130,14 @@ namespace Opm
int j = completion->getJ();
int k = completion->getK();
const int* cpgdim = grid.cartdims;
const int* cpgdim = UgGridHelpers::cartDims(grid);
int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
if (cgit == cartesian_to_compressed.end()) {
OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
<< k << " not found in grid (well = " << well->name() << ')');
if (checkCellExistence_)
OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
<< k << " not found in grid (well = " << well->name() << ')');
continue;
}
int cell = cgit->second;
PerfData pd;
@ -1137,7 +1151,8 @@ namespace Opm
OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
}
std::array<double, 3> cubical = getCubeDim(grid, cell);
const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell];
int dimensions=UgGridHelpers::dimensions(grid);
const double* cell_perm = &permeability[dimensions*dimensions*cell];
pd.well_index = computeWellIndex(radius, cubical, cell_perm, completion->getDiameter());
}
wellperf_data[well_index].push_back(pd);
@ -1151,7 +1166,7 @@ namespace Opm
const int num_wells = well_data.size();
int num_perfs = 0;
assert(grid.dimensions == 3);
assert(=UgGridHelpers::dimensions(grid) == 3);
for (int w = 0; w < num_wells; ++w) {
num_perfs += wellperf_data[w].size();
if (well_data[w].reference_bhp_depth < 0.0) {
@ -1159,7 +1174,8 @@ namespace Opm
double min_depth = 1e100;
int num_wperfs = wellperf_data[w].size();
for (int perf = 0; perf < num_wperfs; ++perf) {
double depth = grid.cell_centroids[3*wellperf_data[w][perf].cell + 2];
double depth = UgGridHelpers::
cellCentroidCoordinate(grid,wellperf_data[w][perf].cell, 2);
min_depth = std::min(min_depth, depth);
}
well_data[w].reference_bhp_depth = min_depth;

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

@ -65,7 +65,7 @@ namespace Opm
/// manage control switching does not exist.
///
/// @param[in] W Existing wells object.
WellsManager(struct Wells* W);
WellsManager(struct Wells* W, bool checkCellExistence=true);
/// Construct from input deck and grid.
/// The permeability argument may be zero if the input contain
@ -73,14 +73,16 @@ namespace Opm
/// order to approximate these by the Peaceman formula.
WellsManager(const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,
const double* permeability);
const double* permeability,
bool checkCellExistence=true);
WellsManager(const Opm::EclipseStateConstPtr eclipseState,
const size_t timeStep,
const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,
const double* permeability);
const double* permeability,
bool checkCellExistence=true);
/// Destructor.
~WellsManager();
@ -135,9 +137,9 @@ namespace Opm
private:
// Disable copying and assignment.
WellsManager(const WellsManager& other);
WellsManager(const WellsManager& other, bool checkCellExistence=true);
WellsManager& operator=(const WellsManager& other);
static void setupCompressedToCartesian(const UnstructuredGrid& grid, std::map<int,int>& cartesian_to_compressed );
static void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed );
void setupWellControls(std::vector<WellConstPtr>& wells, size_t timeStep,
std::vector<std::string>& well_names, const PhaseUsage& phaseUsage);
@ -155,6 +157,7 @@ namespace Opm
// Data
Wells* w_;
WellCollection well_collection_;
bool checkCellExistence_;
};
} // namespace Opm