OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests 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 c1986572c4
commit 8da99e11c5
21 changed files with 1073 additions and 293 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists_files.cmake
View File

@@ -29,6 +29,7 @@
# originally generated with the command:
# find opm -name '*.c*' -printf '\t%p\n' | sort
list (APPEND MAIN_SOURCE_FILES
opm/core/grid/GridHelpers.cpp
opm/core/grid/GridManager.cpp
opm/core/grid/grid.c
opm/core/grid/cart_grid.c
@@ -231,6 +232,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/grid/CellQuadrature.hpp
opm/core/grid/ColumnExtract.hpp
opm/core/grid/FaceQuadrature.hpp
opm/core/grid/GridHelpers.hpp
opm/core/grid/GridManager.hpp
opm/core/grid/cart_grid.h
opm/core/grid/cornerpoint_grid.h
@@ -287,6 +289,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/props/BlackoilPhases.hpp
opm/core/props/BlackoilPropertiesBasic.hpp
opm/core/props/BlackoilPropertiesFromDeck.hpp
opm/core/props/BlackoilPropertiesFromDeck_impl.hpp
opm/core/props/BlackoilPropertiesInterface.hpp
opm/core/props/IncompPropertiesBasic.hpp
opm/core/props/IncompPropertiesFromDeck.hpp
@@ -374,6 +377,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/utility/have_boost_redef.hpp
opm/core/utility/linearInterpolation.hpp
opm/core/utility/miscUtilities.hpp
opm/core/utility/miscUtilities_impl.hpp
opm/core/utility/miscUtilitiesBlackoil.hpp
opm/core/utility/parameters/Parameter.hpp
opm/core/utility/parameters/ParameterGroup.hpp

66
opm/core/grid/GridHelpers.cpp Normal file
View File

@@ -0,0 +1,66 @@
#include "config.h"
#include <opm/core/grid/GridHelpers.hpp>
namespace Opm
{
namespace UgGridHelpers
{
int numCells(const UnstructuredGrid& grid)
{
return grid.number_of_cells;
}
int numFaces(const UnstructuredGrid& grid)
{
return grid.number_of_faces;
}
int dimensions(const UnstructuredGrid& grid)
{
return grid.dimensions;
}
const int* globalCell(const UnstructuredGrid& grid)
{
return grid.global_cell;
}
const int* cartDims(const UnstructuredGrid& grid)
{
return grid.cartdims;
}
const double* beginCellCentroids(const UnstructuredGrid& grid)
{
return grid.cell_centroids;
}
double cellCentroidCoordinate(const UnstructuredGrid& grid, int cell_index,
int coordinate)
{
return grid.cell_centroids[grid.dimensions*cell_index+coordinate];
}
const double* beginFaceCentroids(const UnstructuredGrid& grid)
{
return grid.face_centroids;
}
const double* faceCentroid(const UnstructuredGrid& grid, int face_index)
{
return grid.face_centroids+face_index*grid.dimensions;
}
SparseTableView cell2Faces(const UnstructuredGrid& grid)
{
return SparseTableView(grid.cell_faces, grid.cell_facepos, numCells(grid));
}
double cellVolume(const UnstructuredGrid& grid, int cell_index)
{
return grid.cell_volumes[cell_index];
}
FaceCellTraits<UnstructuredGrid>::Type faceCells(const UnstructuredGrid& grid)
{
return FaceCellsProxy(grid);
}
}
}

223
opm/core/grid/GridHelpers.hpp Normal file
View File

@@ -0,0 +1,223 @@
/*
Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services
Copyright 2014 Statoil AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_CORE_GRIDHELPERS_HEADER_INCLUDED
#define OPM_CORE_GRIDHELPERS_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <boost/range/iterator_range.hpp>
namespace Opm
{
namespace UgGridHelpers
{
/// \brief Allows viewing a sparse table consisting out of C-array
///
/// This class can be used to convert two int array (like they are
/// in UnstructuredGrid for representing the cell to faces mapping
/// as a sparse table object.
class SparseTableView
{
public:
class IntRange : public boost::iterator_range<const int*>
{
public:
typedef boost::iterator_range<const int*> BaseRowType;
typedef BaseRowType::size_type size_type;
typedef int value_type;
IntRange(const int* start, const int* end)
: BaseRowType(start, end)
{}
};
/// \brief The type of the roww.
typedef boost::iterator_range<const int*> row_type;
/// \brief Creates a sparse table view
/// \param data The array with data of the table.
/// \param offset The offsets of the rows. Row i starts
/// at offset[i] and ends a offset[i+1]
/// \param size The number of entries/rows of the table
SparseTableView(int* data, int *offset, std::size_t size)
: data_(data), offset_(offset), size_(size)
{}
/// \brief Get a row of the the table.
/// \param row The row index.
/// \return The corresponding row.
row_type operator[](std::size_t row) const
{
assert(row<=size());
return row_type(data_ + offset_[row], data_ + offset_[row+1]);
}
/// \brief Get the size of the table.
/// \return the number rows.
std::size_t size() const
{
return size_;
}
/// \brief Get the number of non-zero entries.
std::size_t noEntries() const
{
return offset_[size_];
}
private:
/// \brief The array with data of the table.
const int* data_;
/// \brief offset The offsets of the rows.
///
/// Row i starts at offset[i] and ends a offset[i+1]
const int* offset_;
/// \brief The size, i.e. the number of rows.
std::size_t size_;
};
/// \brief Get the number of cells of a grid.
int numCells(const UnstructuredGrid& grid);
/// \brief Get the number of faces of a grid.
int numFaces(const UnstructuredGrid& grid);
/// \brief Get the dimensions of a grid
int dimensions(const UnstructuredGrid& grid);
/// \brief Get the cartesion dimension of the underlying structured grid.
const int* cartDims(const UnstructuredGrid& grid);
/// \brief Get the local to global index mapping.
///
/// The global index is the index of the active cell
/// in the underlying structured grid.
const int* globalCell(const UnstructuredGrid& grid);
/// \brief Get an iterator over the cell centroids positioned at the first cell.
///
/// The return type needs to be usable with the functions increment, and
/// getCoordinate.
const double* beginCellCentroids(const UnstructuredGrid& grid);
/// \brief Get a coordinate of a specific cell centroid.
/// \brief grid The grid.
/// \brief cell_index The index of the specific cell.
/// \breif coordinate The coordinate index.
double cellCentroidCoordinate(const UnstructuredGrid& grid, int cell_index,
int coordinate);
/// \brief Get an iterator over the face centroids positioned at the first cell.
const double* beginFaceCentroids(const UnstructuredGrid& grid);
/// \brief Get a coordinate of a specific face centroid.
/// \brief grid The grid.
/// \brief face_index The index of the specific face.
/// \breif coordinate The coordinate index.
const double* faceCentroid(const UnstructuredGrid& grid, int face_index);
/// \brief Maps the grid type to the associated type of the cell to faces mapping.
///
/// Provides a type named Type.
/// \tparam T The type of the grid.
template<class T>
struct Cell2FacesTraits
{
};
template<>
struct Cell2FacesTraits<UnstructuredGrid>
{
typedef SparseTableView Type;
};
/// \brief Get the cell to faces mapping of a grid.
Cell2FacesTraits<UnstructuredGrid>::Type
cell2Faces(const UnstructuredGrid& grid);
class FaceCellsProxy
{
public:
FaceCellsProxy(const UnstructuredGrid& grid)
: face_cells_(grid.face_cells)
{}
int operator()(int cell_index, int local_index)
{
return face_cells_[2*cell_index+local_index];
}
private:
const int* face_cells_;
};
template<class T>
struct FaceCellTraits
{};
template<>
struct FaceCellTraits<UnstructuredGrid>
{
typedef FaceCellsProxy Type;
};
/// \brief Get the face to cell mapping of a grid.
FaceCellTraits<UnstructuredGrid>::Type faceCells(const UnstructuredGrid& grid);
/// \brief Increment an iterator over an array that reresents a dense row-major
/// matrix with dims columns
/// \param cc The iterator.
/// \param i The nzumber of rows to increment
/// \param dim The number of columns of the matrix.
template<class T>
T* increment(T* cc, int i, int dim)
{
return cc+(i*dim);
}
/// \brief Increment an iterator over an array that reresents a dense row-major
/// matrix with dims columns
/// \param cc The iterator.
/// \param i The nzumber of rows to increment
template<class T>
T increment(const T& t, int i, int)
{
return t+i;
}
/// \brief Get the i-th corrdinate of a centroid.
/// \brief cc The array.
/// \brief i The index of the coordinate.
/// \tparam T The type of the coordinate of the centroid.
template<class T>
double getCoordinate(T* cc, int i)
{
return cc[i];
}
/// \brief Get the i-th corrdinate of an array.
/// \brief cc The array.
/// \brief i The index of the coordinate.
/// \tparam T The type representing the centroid.
/// Has to provide a method center that returns an array with the coordinates.
template<class T>
double getCoordinate(T t, int i)
{
return t->center()[i];
}
} // end namespace UGGridHelpers
} // end namespace OPM
#endif

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

101
opm/core/props/BlackoilPropertiesFromDeck_impl.hpp Normal file
View File

@@ -0,0 +1,101 @@
namespace Opm
{
template<class T>
void BlackoilPropertiesFromDeck::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)
{
if (init_rock){
rock_.init(deck, number_of_cells, global_cell, cart_dims);
}
pvt_.init(deck, 200);
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension, 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() << ").");
}
}
template<class T>
void BlackoilPropertiesFromDeck::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)
{
if(init_rock){
rock_.init(deck, number_of_cells, global_cell, cart_dims);
}
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, number_of_cells, global_cell, begin_cell_centroids, dimension,
sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegUniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
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, number_of_cells, global_cell, begin_cell_centroids, dimension,
sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
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() << ").");
}
}
}

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

15
opm/core/simulator/SimulatorState.cpp
View File

@@ -40,13 +40,18 @@ SimulatorState::vectorApproxEqual(const std::vector<double>& v1,
void
SimulatorState::init(const UnstructuredGrid& g, int num_phases)
{
init(g.number_of_cells, g.number_of_faces, num_phases);
}
void
SimulatorState::init(int number_of_cells, int number_of_faces, int num_phases)
{
num_phases_ = num_phases;
press_.resize(g.number_of_cells, 0.0);
fpress_.resize(g.number_of_faces, 0.0);
flux_.resize(g.number_of_faces, 0.0);
sat_.resize(num_phases_ * g.number_of_cells, 0.0);
for (int cell = 0; cell < g.number_of_cells; ++cell) {
press_.resize(number_of_cells, 0.0);
fpress_.resize(number_of_faces, 0.0);
flux_.resize(number_of_faces, 0.0);
sat_.resize(num_phases_ * number_of_cells, 0.0);
for (int cell = 0; cell < number_of_cells; ++cell) {
// Defaulting the second saturation to 1.0.
// This will usually be oil in a water-oil case,
// gas in an oil-gas case.

2
opm/core/simulator/SimulatorState.hpp
View File

@@ -17,6 +17,8 @@ namespace Opm
virtual void init(const UnstructuredGrid& g, int num_phases);
virtual void init(int number_of_cells, int number_of_phases, int num_phases);
enum ExtremalSat { MinSat, MaxSat };
protected:

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