Merge pull request #1737 from akva2/import_ebos

Import ebos from eWoms
This commit is contained in:
Atgeirr Flø Rasmussen 2019-03-01 10:19:09 +01:00 committed by GitHub
commit 9c6fea72ac
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
37 changed files with 17613 additions and 1 deletions

View File

@ -128,6 +128,7 @@ opm_add_test(test_gatherdeferredlogger
)
# the production oriented general-purpose ECL simulator
opm_add_test(flow
ONLY_COMPILE
ALWAYS_ENABLE
@ -150,8 +151,17 @@ add_test(NAME flow__version
set_tests_properties(flow__version PROPERTIES
PASS_REGULAR_EXPRESSION "${${project}_LABEL}")
# the research oriented general-purpose ECL simulator ("ebos" == &ecl
# &black-&oil &simulator)
opm_add_test(ebos
SOURCES ebos/ebos.cc
EXE_NAME ebos
ONLY_COMPILE)
if(OPM_GRID_FOUND AND HAVE_ECL_INPUT AND HAVE_ECL_OUTPUT)
install(TARGETS ebos DESTINATION bin)
endif()
include(OpmBashCompletion)
opm_add_bash_completion(flow)
opm_add_bash_completion(ebos)

View File

@ -50,6 +50,8 @@ list (APPEND MAIN_SOURCE_FILES
# originally generated with the command:
# find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort
list (APPEND TEST_SOURCE_FILES
tests/test_equil.cc
tests/test_ecl_output.cc
tests/test_blackoil_amg.cpp
tests/test_convergencereport.cpp
tests/test_graphcoloring.cpp
@ -75,6 +77,13 @@ if(MPI_FOUND)
endif()
list (APPEND TEST_DATA_FILES
tests/SUMMARY_DECK_NON_CONSTANT_POROSITY.DATA
tests/equil_base.DATA
tests/equil_capillary.DATA
tests/equil_capillary_overlap.DATA
tests/equil_capillary_swatinit.DATA
tests/equil_deadfluids.DATA
tests/equil_pbvd_and_pdvd.DATA
tests/VFPPROD1
tests/VFPPROD2
tests/msw.data

View File

@ -0,0 +1,253 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Ewoms::AluCartesianIndexMapper
*/
#ifndef EWOMS_ALU_CARTESIAN_INDEX_MAPPER_HH
#define EWOMS_ALU_CARTESIAN_INDEX_MAPPER_HH
#include <dune/grid/common/datahandleif.hh>
#include <dune/grid/utility/persistentcontainer.hh>
#include <array>
#include <memory>
#include <except>
#include <vector>
#include <cassert>
namespace Ewoms {
/*!
* \brief Interface class to access the logical Cartesian grid as used in industry
* standard simulator decks.
*/
template <class Grid>
class AluCartesianIndexMapper
{
public:
// data handle for communicating global ids during load balance and communication
template <class GridView>
class GlobalIndexDataHandle : public Dune::CommDataHandleIF<GlobalIndexDataHandle<GridView>, int>
{
// global id
class GlobalCellIndex
{
public:
GlobalCellIndex()
: idx_(-1)
{}
GlobalCellIndex& operator=(const int index)
{
idx_ = index;
return *this;
}
int index() const
{ return idx_; }
private:
int idx_;
};
typedef typename Dune::PersistentContainer<Grid, GlobalCellIndex> GlobalIndexContainer;
public:
// constructor copying cartesian index to persistent container
GlobalIndexDataHandle(const GridView& gridView,
std::vector<int>& cartesianIndex)
: gridView_(gridView)
, globalIndex_(gridView.grid(), 0)
, cartesianIndex_(cartesianIndex)
{
globalIndex_.resize();
initialize();
}
// constructor copying cartesian index to persistent container
GlobalIndexDataHandle(const GlobalIndexDataHandle& other) = delete ;
// destrcutor writing load balanced cartesian index back to vector
~GlobalIndexDataHandle()
{ finalize(); }
bool contains(int dim, int codim) const
{ return codim == 0; }
bool fixedsize(int dim, int codim) const
{ return true; }
//! \brief loop over all internal data handlers and call gather for
//! given entity
template<class MessageBufferImp, class EntityType>
void gather(MessageBufferImp& buff, const EntityType& element) const
{
int globalIdx = globalIndex_[element].index();
buff.write(globalIdx);
}
//! \brief loop over all internal data handlers and call scatter for
//! given entity
template<class MessageBufferImp, class EntityType>
void scatter(MessageBufferImp& buff, const EntityType& element, size_t n)
{
int globalIdx = -1;
buff.read(globalIdx);
if (globalIdx >= 0)
{
globalIndex_.resize();
globalIndex_[element] = globalIdx;
}
}
//! \brief loop over all internal data handlers and return sum of data
//! size of given entity
template<class EntityType>
size_t size(const EntityType& en) const
{ return 1; }
protected:
// initialize persistent container from given vector
void initialize()
{
auto idx = cartesianIndex_.begin();
auto it = gridView_.template begin<0>();
const auto& endIt = gridView_.template end<0>();
for (; it != endIt; ++it, ++idx)
globalIndex_[*it] = *idx;
}
// update vector from given persistent container
void finalize()
{
std::vector<int> newIndex ;
newIndex.reserve(gridView_.indexSet().size(0));
auto it = gridView_.template begin<0>();
const auto& endIt = gridView_.template end<0>();
for (; it != endIt; ++it)
newIndex.push_back(globalIndex_[*it].index()) ;
cartesianIndex_.swap(newIndex);
}
GridView gridView_;
GlobalIndexContainer globalIndex_;
std::vector<int>& cartesianIndex_;
};
public:
/** \brief dimension of the grid */
static const int dimension = Grid::dimension ;
/** \brief constructor taking grid */
AluCartesianIndexMapper(const Grid& grid,
const std::array<int, dimension>& cartDims,
const std::vector<int>& cartesianIndex)
: grid_(grid)
, cartesianDimensions_(cartDims)
, cartesianIndex_(cartesianIndex)
, cartesianSize_(computeCartesianSize())
{}
/** \brief return Cartesian dimensions, i.e. number of cells in each direction */
const std::array<int, dimension>& cartesianDimensions() const
{ return cartesianDimensions_; }
/** \brief return total number of cells in the logical Cartesian grid */
int cartesianSize() const
{ return cartesianSize_; }
/** \brief return number of cells in the active grid */
int compressedSize() const
{ return cartesianIndex_.size(); }
/** \brief return index of the cells in the logical Cartesian grid */
int cartesianIndex(const int compressedElementIndex) const
{
assert(compressedElementIndex < compressedSize());
return cartesianIndex_[compressedElementIndex];
}
/** \brief return index of the cells in the logical Cartesian grid */
int cartesianIndex(const std::array<int, dimension>& coords) const
{
int cartIndex = coords[0];
int factor = cartesianDimensions()[0];
for (int i=1; i < dimension; ++i) {
cartIndex += coords[i] * factor;
factor *= cartesianDimensions()[i];
}
return cartIndex;
}
/** \brief return Cartesian coordinate, i.e. IJK, for a given cell */
void cartesianCoordinate(const int compressedElementIndex, std::array<int, dimension>& coords) const
{
int gc = cartesianIndex(compressedElementIndex);
if (dimension == 3) {
coords[0] = gc % cartesianDimensions()[0];
gc /= cartesianDimensions()[0];
coords[1] = gc % cartesianDimensions()[1];
coords[2] = gc / cartesianDimensions()[1];
}
else if (dimension == 2) {
coords[0] = gc % cartesianDimensions()[0];
coords[1] = gc / cartesianDimensions()[0];
}
else if (dimension == 1)
coords[0] = gc ;
else
throw std::invalid_argument("cartesianCoordinate not implemented for dimension " + std::to_string(dimension));
}
template <class GridView>
std::unique_ptr<GlobalIndexDataHandle<GridView> > dataHandle(const GridView& gridView)
{
typedef GlobalIndexDataHandle<GridView> DataHandle ;
assert(&grid_ == &gridView.grid());
return std::unique_ptr<DataHandle>(new DataHandle(gridView, cartesianIndex_));
}
protected:
int computeCartesianSize() const
{
int size = cartesianDimensions()[0];
for (int d=1; d<dimension; ++d)
size *= cartesianDimensions()[d];
return size ;
}
const Grid& grid_;
const std::array<int, dimension> cartesianDimensions_;
std::vector<int> cartesianIndex_;
const int cartesianSize_ ;
};
} // end namespace Ewoms
#endif

593
ebos/collecttoiorank.hh Normal file
View File

@ -0,0 +1,593 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef EWOMS_COLLECT_TO_IO_RANK_HH
#define EWOMS_COLLECT_TO_IO_RANK_HH
#include <opm/output/data/Cells.hpp>
#include <opm/output/data/Solution.hpp>
#include <opm/output/data/Wells.hpp>
#include <opm/grid/common/p2pcommunicator.hh>
#include <dune/grid/utility/persistentcontainer.hh>
#include <dune/grid/common/gridenums.hh>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <opm/material/common/Unused.hpp>
#include <dune/grid/common/mcmgmapper.hh>
#include <stdexcept>
namespace Ewoms {
template <class Vanguard>
class CollectDataToIORank
{
public:
typedef typename Vanguard::Grid Grid;
typedef typename Grid::CollectiveCommunication CollectiveCommunication;
// global id
class GlobalCellIndex
{
int globalId_;
int localIndex_;
bool isInterior_;
public:
GlobalCellIndex()
: globalId_(-1)
, localIndex_(-1)
, isInterior_(true)
{}
void setGhost()
{ isInterior_ = false; }
void setId(int globalId)
{ globalId_ = globalId; }
void setIndex(int localIndex)
{ localIndex_ = localIndex; }
int localIndex () const
{ return localIndex_; }
int id () const
{ return globalId_; }
bool isInterior() const
{ return isInterior_; }
};
typedef typename Dune::PersistentContainer<Grid, GlobalCellIndex> GlobalIndexContainer;
static const int dimension = Grid::dimension;
typedef typename Grid::LeafGridView GridView;
typedef GridView AllGridView;
typedef Dune::Point2PointCommunicator<Dune::SimpleMessageBuffer> P2PCommunicatorType;
typedef typename P2PCommunicatorType::MessageBufferType MessageBufferType;
typedef std::vector<GlobalCellIndex> LocalIndexMapType;
typedef std::vector<int> IndexMapType;
typedef std::vector<IndexMapType> IndexMapStorageType;
class DistributeIndexMapping : public P2PCommunicatorType::DataHandleInterface
{
protected:
const std::vector<int>& distributedGlobalIndex_;
IndexMapType& localIndexMap_;
IndexMapStorageType& indexMaps_;
std::map<int, int> globalPosition_;
std::vector<int>& ranks_;
public:
DistributeIndexMapping(const std::vector<int>& globalIndex,
const std::vector<int>& distributedGlobalIndex,
IndexMapType& localIndexMap,
IndexMapStorageType& indexMaps,
std::vector<int>& ranks)
: distributedGlobalIndex_(distributedGlobalIndex)
, localIndexMap_(localIndexMap)
, indexMaps_(indexMaps)
, globalPosition_()
, ranks_(ranks)
{
size_t size = globalIndex.size();
// create mapping globalIndex --> localIndex
for (size_t index = 0; index < size; ++index)
globalPosition_.insert(std::make_pair(globalIndex[index], index));
// we need to create a mapping from local to global
if (!indexMaps_.empty()) {
ranks_.resize(size, -1);
// for the ioRank create a localIndex to index in global state map
IndexMapType& indexMap = indexMaps_.back();
size_t localSize = localIndexMap_.size();
indexMap.resize(localSize);
for (size_t i=0; i<localSize; ++i)
{
int id = distributedGlobalIndex_[localIndexMap_[i]];
indexMap[i] = globalPosition_[id];
ranks_[indexMap[i]] = ioRank;
}
}
}
void pack(int link, MessageBufferType& buffer)
{
// we should only get one link
if (link != 0)
throw std::logic_error("link in method pack is not 0 as execpted");
// pack all interior global cell id's
int size = localIndexMap_.size();
buffer.write(size);
for (int index = 0; index < size; ++index) {
int globalIdx = distributedGlobalIndex_[localIndexMap_[index]];
buffer.write(globalIdx);
}
}
void unpack(int link, MessageBufferType& buffer)
{
// get index map for current link
IndexMapType& indexMap = indexMaps_[link];
assert(!globalPosition_.empty());
// unpack all interior global cell id's
int numCells = 0;
buffer.read(numCells);
indexMap.resize(numCells);
for (int index = 0; index < numCells; ++index) {
int globalId = -1;
buffer.read(globalId);
assert(globalPosition_.find(globalId) != globalPosition_.end());
indexMap[index] = globalPosition_[globalId];
ranks_[indexMap[index]] = link + 1;
}
}
};
enum { ioRank = 0 };
static const bool needsReordering =
!std::is_same<typename Vanguard::Grid, typename Vanguard::EquilGrid>::value;
CollectDataToIORank(const Vanguard& vanguard)
: toIORankComm_()
{
// index maps only have to be build when reordering is needed
if (!needsReordering && !isParallel())
return;
const CollectiveCommunication& comm = vanguard.grid().comm();
{
std::set<int> send, recv;
typedef typename Vanguard::EquilGrid::LeafGridView EquilGridView;
const EquilGridView equilGridView = vanguard.equilGrid().leafGridView();
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
typedef Dune::MultipleCodimMultipleGeomTypeMapper<EquilGridView> EquilElementMapper;
EquilElementMapper equilElemMapper(equilGridView, Dune::mcmgElementLayout());
#else
typedef Dune::MultipleCodimMultipleGeomTypeMapper<EquilGridView, Dune::MCMGElementLayout> EquilElementMapper;
EquilElementMapper equilElemMapper(equilGridView);
#endif
// We need a mapping from local to global grid, here we
// use equilGrid which represents a view on the global grid
const size_t globalSize = vanguard.equilGrid().leafGridView().size(0);
// reserve memory
globalCartesianIndex_.resize(globalSize, -1);
// loop over all elements (global grid) and store Cartesian index
auto elemIt = vanguard.equilGrid().leafGridView().template begin<0>();
const auto& elemEndIt = vanguard.equilGrid().leafGridView().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = equilElemMapper.index(*elemIt);
int cartElemIdx = vanguard.equilCartesianIndexMapper().cartesianIndex(elemIdx);
globalCartesianIndex_[elemIdx] = cartElemIdx;
}
// the I/O rank receives from all other ranks
if (isIORank()) {
for (int i = 0; i < comm.size(); ++i) {
if (i != ioRank)
recv.insert(i);
}
}
else // all other simply send to the I/O rank
send.insert(ioRank);
localIndexMap_.clear();
const size_t gridSize = vanguard.grid().size(0);
localIndexMap_.reserve(gridSize);
// store the local Cartesian index
IndexMapType distributedCartesianIndex;
distributedCartesianIndex.resize(gridSize, -1);
typedef typename Vanguard::GridView LocalGridView;
const LocalGridView localGridView = vanguard.gridView();
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
typedef Dune::MultipleCodimMultipleGeomTypeMapper<LocalGridView> ElementMapper;
ElementMapper elemMapper(localGridView, Dune::mcmgElementLayout());
#else
typedef Dune::MultipleCodimMultipleGeomTypeMapper<LocalGridView, Dune::MCMGElementLayout> ElementMapper;
ElementMapper elemMapper(localGridView);
#endif
// A mapping for the whole grid (including the ghosts) is needed for restarts
auto eIt = localGridView.template begin<0>();
const auto& eEndIt = localGridView.template end<0>();
for (; eIt != eEndIt; ++eIt) {
const auto element = *eIt;
int elemIdx = elemMapper.index(element);
distributedCartesianIndex[elemIdx] = vanguard.cartesianIndex(elemIdx);
// only store interior element for collection
//assert(element.partitionType() == Dune::InteriorEntity);
localIndexMap_.push_back(elemIdx);
}
// insert send and recv linkage to communicator
toIORankComm_.insertRequest(send, recv);
// need an index map for each rank
indexMaps_.clear();
indexMaps_.resize(comm.size());
// distribute global id's to io rank for later association of dof's
DistributeIndexMapping distIndexMapping(globalCartesianIndex_,
distributedCartesianIndex,
localIndexMap_,
indexMaps_,
globalRanks_);
toIORankComm_.exchange(distIndexMapping);
}
}
class PackUnPackCellData : public P2PCommunicatorType::DataHandleInterface
{
const Opm::data::Solution& localCellData_;
Opm::data::Solution& globalCellData_;
const IndexMapType& localIndexMap_;
const IndexMapStorageType& indexMaps_;
public:
PackUnPackCellData(const Opm::data::Solution& localCellData,
Opm::data::Solution& globalCellData,
const IndexMapType& localIndexMap,
const IndexMapStorageType& indexMaps,
size_t globalSize,
bool isIORank)
: localCellData_(localCellData)
, globalCellData_(globalCellData)
, localIndexMap_(localIndexMap)
, indexMaps_(indexMaps)
{
if (isIORank) {
// add missing data to global cell data
for (const auto& pair : localCellData_) {
const std::string& key = pair.first;
std::size_t containerSize = globalSize;
auto OPM_OPTIM_UNUSED ret = globalCellData_.insert(key, pair.second.dim,
std::vector<double>(containerSize),
pair.second.target);
assert(ret.second);
}
MessageBufferType buffer;
pack(0, buffer);
// the last index map is the local one
doUnpack(indexMaps.back(), buffer);
}
}
// pack all data associated with link
void pack(int link, MessageBufferType& buffer)
{
// we should only get one link
if (link != 0)
throw std::logic_error("link in method pack is not 0 as expected");
// write all cell data registered in local state
for (const auto& pair : localCellData_) {
const auto& data = pair.second.data;
// write all data from local data to buffer
write(buffer, localIndexMap_, data);
}
}
void doUnpack(const IndexMapType& indexMap, MessageBufferType& buffer)
{
// we loop over the data as
// its order governs the order the data got received.
for (auto& pair : localCellData_) {
const std::string& key = pair.first;
auto& data = globalCellData_.data(key);
//write all data from local cell data to buffer
read(buffer, indexMap, data);
}
}
// unpack all data associated with link
void unpack(int link, MessageBufferType& buffer)
{ doUnpack(indexMaps_[link], buffer); }
protected:
template <class Vector>
void write(MessageBufferType& buffer,
const IndexMapType& localIndexMap,
const Vector& vector,
unsigned int offset = 0,
unsigned int stride = 1) const
{
unsigned int size = localIndexMap.size();
buffer.write(size);
assert(vector.size() >= stride * size);
for (unsigned int i=0; i<size; ++i)
{
unsigned int index = localIndexMap[i] * stride + offset;
assert(index < vector.size());
buffer.write(vector[index]);
}
}
template <class Vector>
void read(MessageBufferType& buffer,
const IndexMapType& indexMap,
Vector& vector,
unsigned int offset = 0,
unsigned int stride = 1) const
{
unsigned int size = 0;
buffer.read(size);
assert(size == indexMap.size());
for (unsigned int i=0; i<size; ++i) {
unsigned int index = indexMap[i] * stride + offset;
assert(index < vector.size());
buffer.read(vector[index]);
}
}
};
class PackUnPackWellData : public P2PCommunicatorType::DataHandleInterface
{
const Opm::data::Wells& localWellData_;
Opm::data::Wells& globalWellData_;
public:
PackUnPackWellData(const Opm::data::Wells& localWellData,
Opm::data::Wells& globalWellData,
bool isIORank)
: localWellData_(localWellData)
, globalWellData_(globalWellData)
{
if (isIORank) {
MessageBufferType buffer;
pack(0, buffer);
// pass a dummy_link to satisfy virtual class
int dummyLink = -1;
unpack(dummyLink, buffer);
}
}
// pack all data associated with link
void pack(int link, MessageBufferType& buffer)
{
// we should only get one link
if (link != 0)
throw std::logic_error("link in method pack is not 0 as expected");
localWellData_.write(buffer);
}
// unpack all data associated with link
void unpack(int /*link*/, MessageBufferType& buffer)
{ globalWellData_.read(buffer); }
};
class PackUnPackBlockData : public P2PCommunicatorType::DataHandleInterface
{
const std::map<std::pair<std::string, int>, double>& localBlockData_;
std::map<std::pair<std::string, int>, double>& globalBlockValues_;
public:
PackUnPackBlockData(const std::map<std::pair<std::string, int>, double>& localBlockData,
std::map<std::pair<std::string, int>, double>& globalBlockValues,
bool isIORank)
: localBlockData_(localBlockData)
, globalBlockValues_(globalBlockValues)
{
if (isIORank) {
MessageBufferType buffer;
pack(0, buffer);
// pass a dummyLink to satisfy virtual class
int dummyLink = -1;
unpack(dummyLink, buffer);
}
}
// pack all data associated with link
void pack(int link, MessageBufferType& buffer)
{
// we should only get one link
if (link != 0)
throw std::logic_error("link in method pack is not 0 as expected");
// write all block data
unsigned int size = localBlockData_.size();
buffer.write(size);
for (const auto& map : localBlockData_) {
buffer.write(map.first.first);
buffer.write(map.first.second);
buffer.write(map.second);
}
}
// unpack all data associated with link
void unpack(int /*link*/, MessageBufferType& buffer)
{
// read all block data
unsigned int size = 0;
buffer.read(size);
for (size_t i = 0; i < size; ++i) {
std::string name;
int idx;
double data;
buffer.read(name);
buffer.read(idx);
buffer.read(data);
globalBlockValues_[std::make_pair(name, idx)] = data;
}
}
};
// gather solution to rank 0 for EclipseWriter
void collect(const Opm::data::Solution& localCellData,
const std::map<std::pair<std::string, int>, double>& localBlockData,
const Opm::data::Wells& localWellData)
{
globalCellData_ = {};
globalBlockData_.clear();
globalWellData_.clear();
// index maps only have to be build when reordering is needed
if(!needsReordering && !isParallel())
return;
// this also packs and unpacks the local buffers one ioRank
PackUnPackCellData
packUnpackCellData(localCellData,
globalCellData_,
localIndexMap_,
indexMaps_,
numCells(),
isIORank());
if (!isParallel())
// no need to collect anything.
return;
PackUnPackWellData
packUnpackWellData(localWellData,
globalWellData_,
isIORank());
PackUnPackBlockData
packUnpackBlockData(localBlockData,
globalBlockData_,
isIORank());
toIORankComm_.exchange(packUnpackCellData);
toIORankComm_.exchange(packUnpackWellData);
toIORankComm_.exchange(packUnpackBlockData);
#ifndef NDEBUG
// mkae sure every process is on the same page
toIORankComm_.barrier();
#endif
}
const std::map<std::pair<std::string, int>, double>& globalBlockData() const
{ return globalBlockData_; }
const Opm::data::Solution& globalCellData() const
{ return globalCellData_; }
const Opm::data::Wells& globalWellData() const
{ return globalWellData_; }
bool isIORank() const
{ return toIORankComm_.rank() == ioRank; }
bool isParallel() const
{ return toIORankComm_.size() > 1; }
int localIdxToGlobalIdx(unsigned localIdx) const
{
if (!isParallel())
return localIdx;
// the last indexMap is the local one
const IndexMapType& indexMap = indexMaps_.back();
if (indexMap.empty())
throw std::logic_error("index map is not created on this rank");
if (localIdx > indexMap.size())
throw std::logic_error("local index is outside map range");
return indexMap[localIdx];
}
size_t numCells () const
{ return globalCartesianIndex_.size(); }
const std::vector<int>& globalRanks() const
{ return globalRanks_; }
bool isGlobalIdxOnThisRank(unsigned globalIdx) const
{
if (!isParallel())
return true;
// the last indexMap is the local one
const IndexMapType& indexMap = indexMaps_.back();
if (indexMap.empty())
throw std::logic_error("index map is not created on this rank");
return std::find(indexMap.begin(), indexMap.end(), globalIdx) != indexMap.end();
}
protected:
P2PCommunicatorType toIORankComm_;
IndexMapType globalCartesianIndex_;
IndexMapType localIndexMap_;
IndexMapStorageType indexMaps_;
std::vector<int> globalRanks_;
Opm::data::Solution globalCellData_;
std::map<std::pair<std::string, int>, double> globalBlockData_;
Opm::data::Wells globalWellData_;
};
} // end namespace Ewoms
#endif

50
ebos/ebos.cc Normal file
View File

@ -0,0 +1,50 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \brief A general-purpose simulator for ECL decks using the black-oil model.
*/
#include "config.h"
#include <opm/material/common/quad.hpp>
#include <ewoms/common/start.hh>
#include "eclproblem.hh"
BEGIN_PROPERTIES
NEW_TYPE_TAG(EclProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
// Enable experimental features for ebos: ebos is the research simulator of the OPM
// project. If you're looking for a more stable "production quality" simulator, consider
// using `flow`
SET_BOOL_PROP(EclProblem, EnableExperiments, true);
END_PROPERTIES
int main(int argc, char **argv)
{
typedef TTAG(EclProblem) ProblemTypeTag;
return Ewoms::start<ProblemTypeTag>(argc, argv);
}

220
ebos/eclalugridvanguard.hh Normal file
View File

@ -0,0 +1,220 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
* \copydoc Ewoms::EclAluGridVanguard
*/
#ifndef EWOMS_ECL_ALU_GRID_VANGUARD_HH
#define EWOMS_ECL_ALU_GRID_VANGUARD_HH
#include "eclbasevanguard.hh"
#include "alucartesianindexmapper.hh"
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/common/fromtogridfactory.hh>
#include <opm/grid/CpGrid.hpp>
namespace Ewoms {
template <class TypeTag>
class EclAluGridVanguard;
} // namespace Ewoms
BEGIN_PROPERTIES
NEW_TYPE_TAG(EclAluGridVanguard, INHERITS_FROM(EclBaseVanguard));
// declare the properties
SET_TYPE_PROP(EclAluGridVanguard, Vanguard, Ewoms::EclAluGridVanguard<TypeTag>);
SET_TYPE_PROP(EclAluGridVanguard, Grid, Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming>);
SET_TYPE_PROP(EclAluGridVanguard, EquilGrid, Dune::CpGrid);
END_PROPERTIES
namespace Ewoms {
/*!
* \ingroup EclBlackOilSimulator
*
* \brief Helper class for grid instantiation of ECL file-format using problems.
*
* This class uses Dune::ALUGrid as the simulation grid.
*/
template <class TypeTag>
class EclAluGridVanguard : public EclBaseVanguard<TypeTag>
{
friend class EclBaseVanguard<TypeTag>;
typedef EclBaseVanguard<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
public:
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, EquilGrid) EquilGrid;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
private:
typedef Ewoms::AluCartesianIndexMapper<Grid> CartesianIndexMapper;
typedef Dune::CartesianIndexMapper<EquilGrid> EquilCartesianIndexMapper;
static const int dimension = Grid::dimension;
public:
/*!
* \brief Inherit the constructors from the base class.
*/
using EclBaseVanguard<TypeTag>::EclBaseVanguard;
~EclAluGridVanguard()
{
delete cartesianIndexMapper_;
delete equilCartesianIndexMapper_;
delete grid_;
delete equilGrid_;
}
/*!
* \brief Return a reference to the simulation grid.
*/
Grid& grid()
{ return *grid_; }
/*!
* \brief Return a reference to the simulation grid.
*/
const Grid& grid() const
{ return *grid_; }
/*!
* \brief Returns a refefence to the grid which should be used by the EQUIL
* initialization code.
*
* The EQUIL keyword is used to specify the initial condition of the reservoir in
* hydrostatic equilibrium. Since the code which does this is not accepting arbitrary
* DUNE grids (the code is part of the opm-core module), this is not necessarily the
* same as the grid which is used for the actual simulation.
*/
const EquilGrid& equilGrid() const
{ return *equilGrid_; }
/*!
* \brief Indicates that the initial condition has been computed and the memory used
* by the EQUIL grid can be released.
*
* Depending on the implementation, subsequent accesses to the EQUIL grid lead to
* crashes.
*/
void releaseEquilGrid()
{
delete equilCartesianIndexMapper_;
equilCartesianIndexMapper_ = 0;
delete equilGrid_;
equilGrid_ = 0;
}
/*!
* \brief Distribute the simulation grid over multiple processes
*
* (For parallel simulation runs.)
*/
void loadBalance()
{
auto gridView = grid().leafGridView();
auto dataHandle = cartesianIndexMapper_->dataHandle(gridView);
grid().loadBalance(*dataHandle);
// communicate non-interior cells values
grid().communicate(*dataHandle,
Dune::InteriorBorder_All_Interface,
Dune::ForwardCommunication );
}
/*!
* \brief Returns the object which maps a global element index of the simulation grid
* to the corresponding element index of the logically Cartesian index.
*/
const CartesianIndexMapper& cartesianIndexMapper() const
{ return *cartesianIndexMapper_; }
/*!
* \brief Returns mapper from compressed to cartesian indices for the EQUIL grid
*/
const EquilCartesianIndexMapper& equilCartesianIndexMapper() const
{ return *equilCartesianIndexMapper_; }
protected:
void createGrids_()
{
const auto& gridProps = this->eclState().get3DProperties();
const std::vector<double>& porv = gridProps.getDoubleGridProperty("PORV").getData();
// we use separate grid objects: one for the calculation of the initial condition
// via EQUIL and one for the actual simulation. The reason is that the EQUIL code
// cannot cope with arbitrary Dune grids and is also allergic to distributed
// grids.
/////
// create the EQUIL grid
/////
equilGrid_ = new EquilGrid();
equilGrid_->processEclipseFormat(this->eclState().getInputGrid(),
/*isPeriodic=*/false,
/*flipNormals=*/false,
/*clipZ=*/false,
porv);
cartesianCellId_ = equilGrid_->globalCell();
for (unsigned i = 0; i < dimension; ++i)
cartesianDimension_[i] = equilGrid_->logicalCartesianSize()[i];
equilCartesianIndexMapper_ = new EquilCartesianIndexMapper(*equilGrid_);
/////
// create the simulation grid
/////
Dune::FromToGridFactory<Grid> factory;
grid_ = factory.convert(*equilGrid_, cartesianCellId_);
cartesianIndexMapper_ =
new CartesianIndexMapper(*grid_, cartesianDimension_, cartesianCellId_);
}
void filterConnections_()
{
// not handling the removal of completions for this type of grid yet.
}
Grid* grid_;
EquilGrid* equilGrid_;
std::vector<int> cartesianCellId_;
std::array<int,dimension> cartesianDimension_;
CartesianIndexMapper* cartesianIndexMapper_;
EquilCartesianIndexMapper* equilCartesianIndexMapper_;
};
} // namespace Ewoms
#endif

142
ebos/eclbaseaquifermodel.hh Normal file
View File

@ -0,0 +1,142 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
* \copydoc Ewoms::EclBaseAquiferModel
*/
#ifndef EWOMS_ECL_BASE_AQUIFER_MODEL_HH
#define EWOMS_ECL_BASE_AQUIFER_MODEL_HH
#include <ewoms/common/propertysystem.hh>
BEGIN_PROPERTIES
NEW_PROP_TAG(Simulator);
NEW_PROP_TAG(RateVector);
END_PROPERTIES
namespace Ewoms {
/*!
* \ingroup EclBaseAquiferModel
*
* \brief The base class which specifies the API of aquifer models.
*
* This class only provides the API for the actual aquifer model, it does not do
* anything on its own.
*/
template <class TypeTag>
class EclBaseAquiferModel
{
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
public:
EclBaseAquiferModel(Simulator& simulator)
: simulator_(simulator)
{}
/*!
* \brief Called once the problem has been fully initialized and the initial
* condition has been applied.
*/
void initialSolutionApplied()
{ }
/*!
* \brief This method is called when a new episode (report step) starts.
*/
void beginEpisode()
{ }
/*!
* \brief This method is called when a new time step (substep) starts.
*/
void beginTimeStep()
{ }
/*!
* \brief This method is called before each Newton-Raphson iteration.
*/
void beginIteration()
{ }
/*!
* \brief Add the water which enters or leaves the reservoir due to aquifiers.
*/
template <class Context>
void addToSource(RateVector& rate OPM_UNUSED,
const Context& context OPM_UNUSED,
unsigned spaceIdx OPM_UNUSED,
unsigned timeIdx OPM_UNUSED) const
{ }
/*!
* \brief This method is called after each Newton-Raphson successful iteration.
*
* I.e., no exceptions were thrown during the linearization and linear solution
* procedures.
*/
void endIteration()
{ }
/*!
* \brief This method is called after each successful time step (substep).
*
* I.e., all iterations of the time step were successful and the Newton-Raphson
* algorithm converged.
*/
void endTimeStep()
{ }
/*!
* \brief This method is called once an episode (report step) has been finished
* successfully.
*/
void endEpisode()
{ }
/*!
* \brief Write the internal state of the aquifer model to disk using an ad-hoc file
* format.
*/
template <class Restarter>
void serialize(Restarter& res OPM_UNUSED)
{ }
/*!
* \brief Load the internal state of the aquifer model to disk using an ad-hoc file
* format.
*/
template <class Restarter>
void deserialize(Restarter& res OPM_UNUSED)
{ }
protected:
Simulator& simulator_;
};
} // namespace Ewoms
#endif

522
ebos/eclbasevanguard.hh Normal file
View File

@ -0,0 +1,522 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
* \copydoc Ewoms::EclBaseVanguard
*/
#ifndef EWOMS_ECL_BASE_VANGUARD_HH
#define EWOMS_ECL_BASE_VANGUARD_HH
#include <ewoms/io/basevanguard.hh>
#include <ewoms/common/propertysystem.hh>
#include <ewoms/common/parametersystem.hh>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/Parser/ErrorGuard.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
#if HAVE_MPI
#include <mpi.h>
#endif // HAVE_MPI
#include <vector>
#include <unordered_set>
#include <array>
namespace Ewoms {
template <class TypeTag>
class EclBaseVanguard;
}
BEGIN_PROPERTIES
NEW_TYPE_TAG(EclBaseVanguard);
// declare the properties required by the for the ecl simulator vanguard
NEW_PROP_TAG(Grid);
NEW_PROP_TAG(EquilGrid);
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(EclDeckFileName);
NEW_PROP_TAG(OutputDir);
NEW_PROP_TAG(EnableOpmRstFile);
NEW_PROP_TAG(EclStrictParsing);
NEW_PROP_TAG(EclOutputInterval);
NEW_PROP_TAG(IgnoreKeywords);
SET_STRING_PROP(EclBaseVanguard, IgnoreKeywords, "");
SET_STRING_PROP(EclBaseVanguard, EclDeckFileName, "");
SET_INT_PROP(EclBaseVanguard, EclOutputInterval, -1); // use the deck-provided value
SET_BOOL_PROP(EclBaseVanguard, EnableOpmRstFile, false);
SET_BOOL_PROP(EclBaseVanguard, EclStrictParsing, false);
END_PROPERTIES
namespace Ewoms {
/*!
* \ingroup EclBlackOilSimulator
*
* \brief Helper class for grid instantiation of ECL file-format using problems.
*/
template <class TypeTag>
class EclBaseVanguard : public BaseVanguard<TypeTag>
{
typedef BaseVanguard<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
public:
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
protected:
static const int dimension = Grid::dimension;
public:
/*!
* \brief Register the common run-time parameters for all ECL simulator vanguards.
*/
static void registerParameters()
{
EWOMS_REGISTER_PARAM(TypeTag, std::string, EclDeckFileName,
"The name of the file which contains the ECL deck to be simulated");
EWOMS_REGISTER_PARAM(TypeTag, int, EclOutputInterval,
"The number of report steps that ought to be skipped between two writes of ECL results");
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableOpmRstFile,
"Include OPM-specific keywords in the ECL restart file to enable restart of OPM simulators from these files");
EWOMS_REGISTER_PARAM(TypeTag, std::string, IgnoreKeywords,
"List of Eclipse keywords which should be ignored. As a ':' separated string.");
EWOMS_REGISTER_PARAM(TypeTag, bool, EclStrictParsing,
"Use strict mode for parsing - all errors are collected before the applicaton exists.");
}
/*!
* \brief Returns the canonical path to a deck file.
*
* The input can either be the canonical deck file name or the name of the case
* (i.e., without the .DATA extension)
*/
static boost::filesystem::path canonicalDeckPath(const std::string& caseName)
{
const auto fileExists = [](const boost::filesystem::path& f) -> bool
{
if (!boost::filesystem::exists(f))
return false;
if (boost::filesystem::is_regular_file(f))
return true;
return boost::filesystem::is_symlink(f) && boost::filesystem::is_regular_file(boost::filesystem::read_symlink(f));
};
auto simcase = boost::filesystem::path(caseName);
if (fileExists(simcase))
return simcase;
for (const auto& ext : { std::string("data"), std::string("DATA") }) {
if (fileExists(simcase.replace_extension(ext)))
return simcase;
}
throw std::invalid_argument("Cannot find input case '"+caseName+"'");
}
/*!
* \brief Set the wall time which was spend externally to set up the external data structures
*
* i.e., the objects specified via the other setExternal*() methods.
*/
static void setExternalSetupTime(Scalar t)
{ externalSetupTime_ = t; }
/*!
* \brief Returns the wall time required to set up the simulator before it was born.
*/
static Scalar externalSetupTime()
{ return externalSetupTime_; }
/*!
* \brief Set the Opm::EclipseState and the Opm::Deck object which ought to be used
* when the simulator vanguard is instantiated.
*
* This is basically an optimization: In cases where the ECL input deck must be
* examined to decide which simulator ought to be used, this avoids having to parse
* the input twice. When this method is used, the caller is responsible for lifetime
* management of these two objects, i.e., they are not allowed to be deleted as long
* as the simulator vanguard object is alive.
*/
static void setExternalDeck(Opm::Deck* deck, Opm::EclipseState* eclState)
{
externalDeck_ = deck;
externalEclState_ = eclState;
}
/*!
* \brief Create the grid for problem data files which use the ECL file format.
*
* This is the file format used by the commercial ECLiPSE simulator. Usually it uses
* a cornerpoint description of the grid.
*/
EclBaseVanguard(Simulator& simulator)
: ParentType(simulator)
{
int myRank = 0;
#if HAVE_MPI
MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
#endif
std::string fileName = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName);
if (fileName == "")
throw std::runtime_error("No input deck file has been specified as a command line argument,"
" or via '--ecl-deck-file-name=CASE.DATA'");
fileName = canonicalDeckPath(fileName).string();
// compute the base name of the input file name
const char directorySeparator = '/';
long int i;
for (i = fileName.size(); i >= 0; -- i)
if (fileName[i] == directorySeparator)
break;
std::string baseName = fileName.substr(i + 1, fileName.size());
// remove the extension from the input file
for (i = baseName.size(); i >= 0; -- i)
if (baseName[i] == '.')
break;
std::string rawCaseName;
if (i < 0)
rawCaseName = baseName;
else
rawCaseName = baseName.substr(0, i);
// transform the result to ALL_UPPERCASE
caseName_ = rawCaseName;
std::transform(caseName_.begin(), caseName_.end(), caseName_.begin(), ::toupper);
typedef std::pair<std::string, Opm::InputError::Action> ParseModePair;
typedef std::vector<ParseModePair> ParseModePairs;
ParseModePairs tmp;
tmp.emplace_back(Opm::ParseContext::PARSE_RANDOM_SLASH, Opm::InputError::IGNORE);
tmp.emplace_back(Opm::ParseContext::PARSE_MISSING_DIMS_KEYWORD, Opm::InputError::WARN);
tmp.emplace_back(Opm::ParseContext::SUMMARY_UNKNOWN_WELL, Opm::InputError::WARN);
tmp.emplace_back(Opm::ParseContext::SUMMARY_UNKNOWN_GROUP, Opm::InputError::WARN);
Opm::ParseContext parseContext(tmp);
Opm::ErrorGuard errorGuard;
const std::string ignoredKeywords = EWOMS_GET_PARAM(TypeTag, std::string, IgnoreKeywords);
if (ignoredKeywords.size() > 0) {
size_t pos, offset = 0;
while (true) {
pos = ignoredKeywords.find(':', offset);
if (pos == std::string::npos) {
parseContext.ignoreKeyword(ignoredKeywords.substr(offset));
break;
}
parseContext.ignoreKeyword(ignoredKeywords.substr(offset, pos - offset));
offset = pos + 1;
}
}
if (EWOMS_GET_PARAM(TypeTag, bool , EclStrictParsing))
parseContext.update(Opm::InputError::DELAYED_EXIT1);
if (!externalDeck_) {
if (myRank == 0)
std::cout << "Reading the deck file '" << fileName << "'" << std::endl;
Opm::Parser parser;
internalDeck_.reset(new Opm::Deck(parser.parseFile(fileName , parseContext, errorGuard)));
internalEclState_.reset(new Opm::EclipseState(*internalDeck_, parseContext, errorGuard));
deck_ = &(*internalDeck_);
eclState_ = &(*internalEclState_);
}
else {
assert(externalDeck_);
assert(externalEclState_);
deck_ = externalDeck_;
eclState_ = externalEclState_;
}
if (!externalEclSchedule_) {
// create the schedule object. Note that if eclState is supposed to represent
// the internalized version of the deck, this constitutes a layering
// violation.
internalEclSchedule_.reset(new Opm::Schedule(*deck_, *eclState_, parseContext, errorGuard));
eclSchedule_ = &(*internalEclSchedule_);
}
else
eclSchedule_ = externalEclSchedule_;
if (!externalEclSummaryConfig_) {
// create the schedule object. Note that if eclState is supposed to represent
// the internalized version of the deck, this constitutes a layering
// violation.
internalEclSummaryConfig_.reset(new Opm::SummaryConfig(*deck_,
*eclSchedule_,
eclState_->getTableManager(),
parseContext,
errorGuard));
eclSummaryConfig_ = &(*internalEclSummaryConfig_);
}
else
eclSummaryConfig_ = externalEclSummaryConfig_;
if (errorGuard) {
errorGuard.dump();
errorGuard.clear();
throw std::runtime_error("Unrecoverable errors were encountered while loading input.");
}
// Possibly override IOConfig setting for how often RESTART files should get
// written to disk (every N report step)
int outputInterval = EWOMS_GET_PARAM(TypeTag, int, EclOutputInterval);
if (outputInterval >= 0)
eclState_->getRestartConfig().overrideRestartWriteInterval(outputInterval);
asImp_().createGrids_();
asImp_().filterConnections_();
asImp_().updateOutputDir_();
asImp_().finalizeInit_();
}
/*!
* \brief Return a reference to the parsed ECL deck.
*/
const Opm::Deck& deck() const
{ return *deck_; }
Opm::Deck& deck()
{ return *deck_; }
/*!
* \brief Return a reference to the internalized ECL deck.
*/
const Opm::EclipseState& eclState() const
{ return *eclState_; }
Opm::EclipseState& eclState()
{ return *eclState_; }
/*!
* \brief Return a reference to the object that managages the ECL schedule.
*/
const Opm::Schedule& schedule() const
{ return *eclSchedule_; }
Opm::Schedule& schedule()
{ return *eclSchedule_; }
/*!
* \brief Set the schedule object.
*
* The lifetime of this object is not managed by the vanguard, i.e., the object must
* stay valid until after the vanguard gets destroyed.
*/
static void setExternalSchedule(Opm::Schedule* schedule)
{ externalEclSchedule_ = schedule; }
/*!
* \brief Return a reference to the object that determines which quantities ought to
* be put into the ECL summary output.
*/
const Opm::SummaryConfig& summaryConfig() const
{ return *eclSummaryConfig_; }
/*!
* \brief Set the summary configuration object.
*
* The lifetime of this object is not managed by the vanguard, i.e., the object must
* stay valid until after the vanguard gets destroyed.
*/
static void setExternalSummaryConfig(Opm::SummaryConfig* summaryConfig)
{ externalEclSummaryConfig_ = summaryConfig; }
/*!
* \brief Returns the name of the case.
*
* i.e., the all-uppercase version of the file name from which the
* deck is loaded with the ".DATA" suffix removed.
*/
const std::string& caseName() const
{ return caseName_; }
/*!
* \brief Returns the number of logically Cartesian cells in each direction
*/
const std::array<int, dimension>& cartesianDimensions() const
{ return asImp_().cartesianIndexMapper().cartesianDimensions(); }
/*!
* \brief Returns the overall number of cells of the logically Cartesian grid
*/
int cartesianSize() const
{ return asImp_().cartesianIndexMapper().cartesianSize(); }
/*!
* \brief Returns the overall number of cells of the logically EquilCartesian grid
*/
int equilCartesianSize() const
{ return asImp_().equilCartesianIndexMapper().cartesianSize(); }
/*!
* \brief Returns the Cartesian cell id for identifaction with ECL data
*/
unsigned cartesianIndex(unsigned compressedCellIdx) const
{ return asImp_().cartesianIndexMapper().cartesianIndex(compressedCellIdx); }
/*!
* \brief Return the index of the cells in the logical Cartesian grid
*/
unsigned cartesianIndex(const std::array<int,dimension>& coords) const
{
unsigned cartIndex = coords[0];
int factor = cartesianDimensions()[0];
for (unsigned i = 1; i < dimension; ++i) {
cartIndex += coords[i]*factor;
factor *= cartesianDimensions()[i];
}
return cartIndex;
}
/*!
* \brief Extract Cartesian index triplet (i,j,k) of an active cell.
*
* \param [in] cellIdx Active cell index.
* \param [out] ijk Cartesian index triplet
*/
void cartesianCoordinate(unsigned cellIdx, std::array<int,3>& ijk) const
{ return asImp_().cartesianIndexMapper().cartesianCoordinate(cellIdx, ijk); }
/*!
* \brief Returns the Cartesian cell id given an element index for the grid used for equilibration
*/
unsigned equilCartesianIndex(unsigned compressedEquilCellIdx) const
{ return asImp_().equilCartesianIndexMapper().cartesianIndex(compressedEquilCellIdx); }
/*!
* \brief Extract Cartesian index triplet (i,j,k) of an active cell of the grid used for EQUIL.
*
* \param [in] cellIdx Active cell index.
* \param [out] ijk Cartesian index triplet
*/
void equilCartesianCoordinate(unsigned cellIdx, std::array<int,3>& ijk) const
{ return asImp_().equilCartesianIndexMapper().cartesianCoordinate(cellIdx, ijk); }
/*!
* \brief Return the names of the wells which do not penetrate any cells on the local
* process.
*
* This is a kludge around the fact that for distributed grids, not all wells are
* seen by all proccesses.
*/
std::unordered_set<std::string> defunctWellNames() const
{ return std::unordered_set<std::string>(); }
private:
void updateOutputDir_()
{
// update the location for output
std::string outputDir = EWOMS_GET_PARAM(TypeTag, std::string, OutputDir);
auto& ioConfig = eclState_->getIOConfig();
if (outputDir == "")
// If no output directory parameter is specified, use the output directory
// which Opm::IOConfig thinks that should be used. Normally this is the
// directory in which the input files are located.
outputDir = ioConfig.getOutputDir();
// ensure that the output directory exists and that it is a directory
if (!boost::filesystem::is_directory(outputDir)) {
try {
boost::filesystem::create_directories(outputDir);
}
catch (...) {
throw std::runtime_error("Creation of output directory '"+outputDir+"' failed\n");
}
}
// specify the directory output. This is not a very nice mechanism because
// the eclState is supposed to be immutable here, IMO.
ioConfig.setOutputDir(outputDir);
ioConfig.setEclCompatibleRST(!EWOMS_GET_PARAM(TypeTag, bool, EnableOpmRstFile));
}
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
const Implementation& asImp_() const
{ return *static_cast<const Implementation*>(this); }
std::string caseName_;
static Scalar externalSetupTime_;
static Opm::Deck* externalDeck_;
static Opm::EclipseState* externalEclState_;
static Opm::Schedule* externalEclSchedule_;
static Opm::SummaryConfig* externalEclSummaryConfig_;
std::unique_ptr<Opm::Deck> internalDeck_;
std::unique_ptr<Opm::EclipseState> internalEclState_;
std::unique_ptr<Opm::Schedule> internalEclSchedule_;
std::unique_ptr<Opm::SummaryConfig> internalEclSummaryConfig_;
// these two attributes point either to the internal or to the external version of the
// Deck and EclipsState objects.
Opm::Deck* deck_;
Opm::EclipseState* eclState_;
Opm::Schedule* eclSchedule_;
Opm::SummaryConfig* eclSummaryConfig_;
};
template <class TypeTag>
typename EclBaseVanguard<TypeTag>::Scalar EclBaseVanguard<TypeTag>::externalSetupTime_ = 0.0;
template <class TypeTag>
Opm::Deck* EclBaseVanguard<TypeTag>::externalDeck_ = nullptr;
template <class TypeTag>
Opm::EclipseState* EclBaseVanguard<TypeTag>::externalEclState_;
template <class TypeTag>
Opm::Schedule* EclBaseVanguard<TypeTag>::externalEclSchedule_ = nullptr;
template <class TypeTag>
Opm::SummaryConfig* EclBaseVanguard<TypeTag>::externalEclSummaryConfig_ = nullptr;
} // namespace Ewoms
#endif

292
ebos/eclcpgridvanguard.hh Normal file
View File

@ -0,0 +1,292 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
* \copydoc Ewoms::EclCpGridVanguard
*/
#ifndef EWOMS_ECL_CP_GRID_VANGUARD_HH
#define EWOMS_ECL_CP_GRID_VANGUARD_HH
#include "eclbasevanguard.hh"
#include "ecltransmissibility.hh"
#include "femcpgridcompat.hh"
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/cpgrid/GridHelpers.hpp>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/common/version.hh>
namespace Ewoms {
template <class TypeTag>
class EclCpGridVanguard;
}
BEGIN_PROPERTIES
NEW_TYPE_TAG(EclCpGridVanguard, INHERITS_FROM(EclBaseVanguard));
// declare the properties
SET_TYPE_PROP(EclCpGridVanguard, Vanguard, Ewoms::EclCpGridVanguard<TypeTag>);
SET_TYPE_PROP(EclCpGridVanguard, Grid, Dune::CpGrid);
SET_TYPE_PROP(EclCpGridVanguard, EquilGrid, typename GET_PROP_TYPE(TypeTag, Grid));
END_PROPERTIES
namespace Ewoms {
/*!
* \ingroup EclBlackOilSimulator
*
* \brief Helper class for grid instantiation of ECL file-format using problems.
*
* This class uses Dune::CpGrid as the simulation grid.
*/
template <class TypeTag>
class EclCpGridVanguard : public EclBaseVanguard<TypeTag>
{
friend class EclBaseVanguard<TypeTag>;
typedef EclBaseVanguard<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
public:
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, EquilGrid) EquilGrid;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
private:
typedef Dune::CartesianIndexMapper<Grid> CartesianIndexMapper;
public:
/*!
* \brief Inherit the constructors from the base class.
*/
using EclBaseVanguard<TypeTag>::EclBaseVanguard;
~EclCpGridVanguard()
{
delete cartesianIndexMapper_;
delete equilCartesianIndexMapper_;
delete grid_;
delete equilGrid_;
delete globalTrans_;
}
/*!
* \brief Return a reference to the simulation grid.
*/
Grid& grid()
{ return *grid_; }
/*!
* \brief Return a reference to the simulation grid.
*/
const Grid& grid() const
{ return *grid_; }
/*!
* \brief Returns a refefence to the grid which should be used by the EQUIL
* initialization code.
*
* The EQUIL keyword is used to specify the initial condition of the reservoir in
* hydrostatic equilibrium. Since the code which does this is not accepting arbitrary
* DUNE grids (the code is part of the opm-core module), this is not necessarily the
* same as the grid which is used for the actual simulation.
*/
const EquilGrid& equilGrid() const
{ return *equilGrid_; }
/*!
* \brief Indicates that the initial condition has been computed and the memory used
* by the EQUIL grid can be released.
*
* Depending on the implementation, subsequent accesses to the EQUIL grid lead to
* crashes.
*/
void releaseEquilGrid()
{
delete equilGrid_;
equilGrid_ = 0;
delete equilCartesianIndexMapper_;
equilCartesianIndexMapper_ = 0;
}
/*!
* \brief Distribute the simulation grid over multiple processes
*
* (For parallel simulation runs.)
*/
void loadBalance()
{
#if HAVE_MPI
int mpiRank = 0;
int mpiSize = 1;
MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank);
MPI_Comm_size(MPI_COMM_WORLD, &mpiSize);
if (mpiSize > 1) {
// the CpGrid's loadBalance() method likes to have the transmissibilities as
// its edge weights. since this is (kind of) a layering violation and
// transmissibilities are relatively expensive to compute, we only do it if
// more than a single process is involved in the simulation.
cartesianIndexMapper_ = new CartesianIndexMapper(*grid_);
globalTrans_ = new EclTransmissibility<TypeTag>(*this);
globalTrans_->update();
// convert to transmissibility for faces
// TODO: grid_->numFaces() is not generic. use grid_->size(1) instead? (might
// not work)
const auto& gridView = grid_->leafGridView();
unsigned numFaces = grid_->numFaces();
std::vector<double> faceTrans(numFaces, 0.0);
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
#else
ElementMapper elemMapper(this->gridView());
#endif
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
const auto& elem = *elemIt;
auto isIt = gridView.ibegin(elem);
const auto& isEndIt = gridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
const auto& is = *isIt;
if (!is.neighbor())
continue;
unsigned I = elemMapper.index(is.inside());
unsigned J = elemMapper.index(is.outside());
// FIXME (?): this is not portable!
unsigned faceIdx = is.id();
faceTrans[faceIdx] = globalTrans_->transmissibility(I, J);
}
}
//distribute the grid and switch to the distributed view.
{
const auto wells = this->schedule().getWells();
defunctWellNames_ = std::get<1>(grid_->loadBalance(&wells, faceTrans.data()));
}
grid_->switchToDistributedView();
delete cartesianIndexMapper_;
cartesianIndexMapper_ = nullptr;
}
#endif
cartesianIndexMapper_ = new CartesianIndexMapper(*grid_);
this->updateGridView_();
}
/*!
* \brief Free the memory occupied by the global transmissibility object.
*
* After writing the initial solution, this array should not be necessary anymore.
*/
void releaseGlobalTransmissibilities()
{
delete globalTrans_;
globalTrans_ = nullptr;
}
/*!
* \brief Returns the object which maps a global element index of the simulation grid
* to the corresponding element index of the logically Cartesian index.
*/
const CartesianIndexMapper& cartesianIndexMapper() const
{ return *cartesianIndexMapper_; }
/*!
* \brief Returns mapper from compressed to cartesian indices for the EQUIL grid
*/
const CartesianIndexMapper& equilCartesianIndexMapper() const
{ return *equilCartesianIndexMapper_; }
std::unordered_set<std::string> defunctWellNames() const
{ return defunctWellNames_; }
const EclTransmissibility<TypeTag>& globalTransmissibility() const
{ return *globalTrans_; }
void releaseGlobalTransmissibility()
{
delete globalTrans_;
globalTrans_ = nullptr;
}
protected:
void createGrids_()
{
const auto& gridProps = this->eclState().get3DProperties();
const std::vector<double>& porv = gridProps.getDoubleGridProperty("PORV").getData();
grid_ = new Dune::CpGrid();
grid_->processEclipseFormat(this->eclState().getInputGrid(),
/*isPeriodic=*/false,
/*flipNormals=*/false,
/*clipZ=*/false,
porv,
this->eclState().getInputNNC());
// we use separate grid objects: one for the calculation of the initial condition
// via EQUIL and one for the actual simulation. The reason is that the EQUIL code
// is allergic to distributed grids and the simulation grid is distributed before
// the initial condition is calculated.
// After loadbalance grid_ will contain a global and distribute view.
// equilGrid_being a shallow copy only the global view.
equilGrid_ = new Dune::CpGrid(*grid_);
equilCartesianIndexMapper_ = new CartesianIndexMapper(*equilGrid_);
globalTrans_ = nullptr;
}
// removing some connection located in inactive grid cells
void filterConnections_()
{
assert(grid_);
Grid grid = *grid_;
grid.switchToGlobalView();
const auto eclipseGrid = Opm::UgGridHelpers::createEclipseGrid(grid, this->eclState().getInputGrid());
this->schedule().filterConnections(eclipseGrid);
}
Grid* grid_;
EquilGrid* equilGrid_;
CartesianIndexMapper* cartesianIndexMapper_;
CartesianIndexMapper* equilCartesianIndexMapper_;
EclTransmissibility<TypeTag>* globalTrans_;
std::unordered_set<std::string> defunctWellNames_;
};
} // namespace Ewoms
#endif

View File

@ -0,0 +1,104 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Ewoms::EclDummyGradientCalculator
*/
#ifndef EWOMS_ECL_DUMMY_GRADIENT_CALCULATOR_HH
#define EWOMS_ECL_DUMMY_GRADIENT_CALCULATOR_HH
#include <ewoms/disc/common/fvbaseproperties.hh>
#include <opm/material/common/Exceptions.hpp>
#include <dune/common/fvector.hh>
namespace Ewoms {
/*!
* \ingroup EclBlackOilSimulator
*
* \brief This is a "dummy" gradient calculator which does not do anything.
*
* The ECL blackoil simulator does not need any gradients: Volume fluxes are calculated
* via pressure differences instead of pressure gradients (i.e., transmissibilities
* instead of permeabilities), and an energy equation and molecular diffusion are not
* supported.
*/
template<class TypeTag>
class EclDummyGradientCalculator
{
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
enum { dimWorld = GridView::dimensionworld };
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
public:
static void registerParameters()
{ }
template <bool prepareValues = true, bool prepareGradients = true>
void prepare(const ElementContext& elemCtx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
{ }
template <class QuantityCallback, class QuantityType = Scalar>
QuantityType calculateValue(const ElementContext& elemCtx OPM_UNUSED,
unsigned fapIdx OPM_UNUSED,
const QuantityCallback& quantityCallback OPM_UNUSED) const
{
throw std::logic_error("Generic values are not supported by the ECL black-oil simulator");
}
template <class QuantityCallback>
void calculateGradient(DimVector& quantityGrad OPM_UNUSED,
const ElementContext& elemCtx OPM_UNUSED,
unsigned fapIdx OPM_UNUSED,
const QuantityCallback& quantityCallback OPM_UNUSED) const
{
throw std::logic_error("Generic gradients are not supported by the ECL black-oil simulator");
}
template <class QuantityCallback>
Scalar calculateBoundaryValue(const ElementContext& elemCtx OPM_UNUSED,
unsigned fapIdx OPM_UNUSED,
const QuantityCallback& quantityCallback OPM_UNUSED)
{
throw std::logic_error("Generic boundary values are not supported by the ECL black-oil simulator");
}
template <class QuantityCallback>
void calculateBoundaryGradient(DimVector& quantityGrad OPM_UNUSED,
const ElementContext& elemCtx OPM_UNUSED,
unsigned fapIdx OPM_UNUSED,
const QuantityCallback& quantityCallback OPM_UNUSED) const
{
throw std::logic_error("Generic boundary gradients are not supported by the ECL black-oil simulator");
}
};
} // namespace Ewoms
#endif

186
ebos/eclequilinitializer.hh Normal file
View File

@ -0,0 +1,186 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/**
* \file
*
* \copydoc Ewoms::EclEquilInitializer
*/
#ifndef EWOMS_ECL_EQUIL_INITIALIZER_HH
#define EWOMS_ECL_EQUIL_INITIALIZER_HH
#include "equil/initstateequil.hh"
#include <ewoms/common/propertysystem.hh>
#include <opm/material/fluidstates/BlackOilFluidState.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <vector>
BEGIN_PROPERTIES
NEW_PROP_TAG(Simulator);
NEW_PROP_TAG(FluidSystem);
NEW_PROP_TAG(GridView);
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(MaterialLaw);
NEW_PROP_TAG(EnableTemperature);
NEW_PROP_TAG(EnableEnergy);
END_PROPERTIES
namespace Ewoms {
/*!
* \ingroup EclBlackOilSimulator
*
* \brief Computes the initial condition based on the EQUIL keyword from ECL.
*
* So far, it uses the "initStateEquil()" function from opm-core. Since this method is
* very much glued into the opm-core data structures, it should be reimplemented in the
* medium to long term for some significant memory savings and less significant
* performance improvements.
*/
template <class TypeTag>
class EclEquilInitializer
{
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum { numPhases = FluidSystem::numPhases };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
enum { numComponents = FluidSystem::numComponents };
enum { oilCompIdx = FluidSystem::oilCompIdx };
enum { gasCompIdx = FluidSystem::gasCompIdx };
enum { waterCompIdx = FluidSystem::waterCompIdx };
enum { dimWorld = GridView::dimensionworld };
enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
public:
// NB: setting the enableEnergy argument to true enables storage of enthalpy and
// internal energy!
typedef Opm::BlackOilFluidState<Scalar,
FluidSystem,
enableTemperature,
enableEnergy,
Indices::gasEnabled,
Indices::numPhases
> ScalarFluidState;
template <class EclMaterialLawManager>
EclEquilInitializer(const Simulator& simulator,
EclMaterialLawManager& materialLawManager)
: simulator_(simulator)
{
const auto& vanguard = simulator.vanguard();
const auto& eclState = vanguard.eclState();
unsigned numElems = vanguard.grid().size(0);
unsigned numCartesianElems = vanguard.cartesianSize();
EQUIL::DeckDependent::InitialStateComputer<TypeTag> initialState(materialLawManager,
eclState,
vanguard.grid(),
simulator.problem().gravity()[dimWorld - 1]);
// copy the result into the array of initial fluid states
initialFluidStates_.resize(numCartesianElems);
for (unsigned int elemIdx = 0; elemIdx < numElems; ++elemIdx) {
unsigned cartesianElemIdx = vanguard.cartesianIndex(elemIdx);
auto& fluidState = initialFluidStates_[cartesianElemIdx];
// get the PVT region index of the current element
unsigned regionIdx = simulator_.problem().pvtRegionIndex(elemIdx);
fluidState.setPvtRegionIndex(regionIdx);
// set the phase saturations
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (FluidSystem::phaseIsActive(phaseIdx))
fluidState.setSaturation(phaseIdx, initialState.saturation()[phaseIdx][elemIdx]);
else if (Indices::numPhases == 3)
fluidState.setSaturation(phaseIdx, 0.0);
}
if (FluidSystem::enableDissolvedGas())
fluidState.setRs(initialState.rs()[elemIdx]);
else if (Indices::gasEnabled)
fluidState.setRs(0.0);
if (FluidSystem::enableVaporizedOil())
fluidState.setRv(initialState.rv()[elemIdx]);
else if (Indices::gasEnabled)
fluidState.setRv(0.0);
// set the temperature.
if (enableTemperature || enableEnergy)
fluidState.setTemperature(initialState.temperature()[elemIdx]);
// set the phase pressures, invB factor and density
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
fluidState.setPressure(phaseIdx, initialState.press()[phaseIdx][elemIdx]);
const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, regionIdx);
fluidState.setInvB(phaseIdx, b);
const auto& rho = FluidSystem::density(fluidState, phaseIdx, regionIdx);
fluidState.setDensity(phaseIdx, rho);
}
}
}
/*!
* \brief Return the initial thermodynamic state which should be used as the initial
* condition.
*
* This is supposed to correspond to hydrostatic conditions.
*/
const ScalarFluidState& initialFluidState(unsigned elemIdx) const
{
const auto& vanguard = simulator_.vanguard();
unsigned cartesianElemIdx = vanguard.cartesianIndex(elemIdx);
return initialFluidStates_[cartesianElemIdx];
}
protected:
const Simulator& simulator_;
std::vector<ScalarFluidState> initialFluidStates_;
};
} // namespace Ewoms
#endif

493
ebos/eclfluxmodule.hh Normal file
View File

@ -0,0 +1,493 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \brief This file contains the flux module which is used for ECL problems
*
* This approach to fluxes is very specific to two-point flux approximation and applies
* what the Eclipse Technical Description calls the "NEWTRAN" transmissibility approach.
*/
#ifndef EWOMS_ECL_FLUX_MODULE_HH
#define EWOMS_ECL_FLUX_MODULE_HH
#include <ewoms/disc/common/fvbaseproperties.hh>
#include <ewoms/models/blackoil/blackoilproperties.hh>
#include <ewoms/common/signum.hh>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
BEGIN_PROPERTIES
NEW_PROP_TAG(MaterialLaw);
END_PROPERTIES
namespace Ewoms {
template <class TypeTag>
class EclTransIntensiveQuantities;
template <class TypeTag>
class EclTransExtensiveQuantities;
template <class TypeTag>
class EclTransBaseProblem;
/*!
* \ingroup EclBlackOilSimulator
* \brief Specifies a flux module which uses ECL transmissibilities.
*/
template <class TypeTag>
struct EclTransFluxModule
{
typedef EclTransIntensiveQuantities<TypeTag> FluxIntensiveQuantities;
typedef EclTransExtensiveQuantities<TypeTag> FluxExtensiveQuantities;
typedef EclTransBaseProblem<TypeTag> FluxBaseProblem;
/*!
* \brief Register all run-time parameters for the flux module.
*/
static void registerParameters()
{ }
};
/*!
* \ingroup EclBlackOilSimulator
* \brief Provides the defaults for the parameters required by the
* transmissibility based volume flux calculation.
*/
template <class TypeTag>
class EclTransBaseProblem
{ };
/*!
* \ingroup EclBlackOilSimulator
* \brief Provides the intensive quantities for the ECL flux module
*/
template <class TypeTag>
class EclTransIntensiveQuantities
{
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
protected:
void update_(const ElementContext& elemCtx OPM_UNUSED, unsigned dofIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
{ }
};
/*!
* \ingroup EclBlackOilSimulator
* \brief Provides the ECL flux module
*/
template <class TypeTag>
class EclTransExtensiveQuantities
{
typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
enum { dimWorld = GridView::dimensionworld };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
enum { numPhases = FluidSystem::numPhases };
enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
typedef Dune::FieldVector<Evaluation, dimWorld> EvalDimVector;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
public:
/*!
* \brief Return the intrinsic permeability tensor at a face [m^2]
*/
const DimMatrix& intrinsicPermeability() const
{
throw std::invalid_argument("The ECL transmissibility module does not provide an explicit intrinsic permeability");
}
/*!
* \brief Return the pressure potential gradient of a fluid phase at the
* face's integration point [Pa/m]
*
* \param phaseIdx The index of the fluid phase
*/
const EvalDimVector& potentialGrad(unsigned phaseIdx OPM_UNUSED) const
{
throw std::invalid_argument("The ECL transmissibility module does not provide explicit potential gradients");
}
/*!
* \brief Return the gravity corrected pressure difference between the interior and
* the exterior of a face.
*
* \param phaseIdx The index of the fluid phase
*/
const Evaluation& pressureDifference(unsigned phaseIdx) const
{ return pressureDifference_[phaseIdx]; }
/*!
* \brief Return the filter velocity of a fluid phase at the face's integration point
* [m/s]
*
* \param phaseIdx The index of the fluid phase
*/
const EvalDimVector& filterVelocity(unsigned phaseIdx OPM_UNUSED) const
{
throw std::invalid_argument("The ECL transmissibility module does not provide explicit filter velocities");
}
/*!
* \brief Return the volume flux of a fluid phase at the face's integration point
* \f$[m^3/s / m^2]\f$
*
* This is the fluid volume of a phase per second and per square meter of face
* area.
*
* \param phaseIdx The index of the fluid phase
*/
const Evaluation& volumeFlux(unsigned phaseIdx) const
{ return volumeFlux_[phaseIdx]; }
protected:
/*!
* \brief Returns the local index of the degree of freedom in which is
* in upstream direction.
*
* i.e., the DOF which exhibits a higher effective pressure for
* the given phase.
*/
unsigned upstreamIndex_(unsigned phaseIdx) const
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
return upIdx_[phaseIdx];
}
/*!
* \brief Returns the local index of the degree of freedom in which is
* in downstream direction.
*
* i.e., the DOF which exhibits a lower effective pressure for the
* given phase.
*/
unsigned downstreamIndex_(unsigned phaseIdx) const
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
return dnIdx_[phaseIdx];
}
void updateSolvent(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
{ asImp_().updateVolumeFluxTrans(elemCtx, scvfIdx, timeIdx); }
void updatePolymer(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
{ asImp_().updateShearMultipliers(elemCtx, scvfIdx, timeIdx); }
/*!
* \brief Update the required gradients for interior faces
*/
void calculateGradients_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
{
Opm::Valgrind::SetUndefined(*this);
const auto& problem = elemCtx.problem();
const auto& stencil = elemCtx.stencil(timeIdx);
const auto& scvf = stencil.interiorFace(scvfIdx);
interiorDofIdx_ = scvf.interiorIndex();
exteriorDofIdx_ = scvf.exteriorIndex();
assert(interiorDofIdx_ != exteriorDofIdx_);
unsigned I = stencil.globalSpaceIndex(interiorDofIdx_);
unsigned J = stencil.globalSpaceIndex(exteriorDofIdx_);
Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx_, exteriorDofIdx_);
Scalar faceArea = scvf.area();
Scalar thpres = problem.thresholdPressure(I, J);
// estimate the gravity correction: for performance reasons we use a simplified
// approach for this flux module that assumes that gravity is constant and always
// acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx_, timeIdx);
const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx_, timeIdx);
// this is quite hacky because the dune grid interface does not provide a
// cellCenterDepth() method (so we ask the problem to provide it). The "good"
// solution would be to take the Z coordinate of the element centroids, but since
// ECL seems to like to be inconsistent on that front, it needs to be done like
// here...
Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx_, timeIdx);
Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx_, timeIdx);
// the distances from the DOF's depths. (i.e., the additional depth of the
// exterior DOF)
Scalar distZ = zIn - zEx;
for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
// check shortcut: if the mobility of the phase is zero in the interior as
// well as the exterior DOF, we can skip looking at the phase.
if (intQuantsIn.mobility(phaseIdx) <= 0.0 &&
intQuantsEx.mobility(phaseIdx) <= 0.0)
{
upIdx_[phaseIdx] = interiorDofIdx_;
dnIdx_[phaseIdx] = exteriorDofIdx_;
pressureDifference_[phaseIdx] = 0.0;
volumeFlux_[phaseIdx] = 0.0;
continue;
}
// do the gravity correction: compute the hydrostatic pressure for the
// external at the depth of the internal one
const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
Evaluation rhoAvg = (rhoIn + rhoEx)/2;
const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx));
pressureExterior += rhoAvg*(distZ*g);
pressureDifference_[phaseIdx] = pressureExterior - pressureInterior;
// decide the upstream index for the phase. for this we make sure that the
// degree of freedom which is regarded upstream if both pressures are equal
// is always the same: if the pressure is equal, the DOF with the lower
// global index is regarded to be the upstream one.
if (pressureDifference_[phaseIdx] > 0.0) {
upIdx_[phaseIdx] = exteriorDofIdx_;
dnIdx_[phaseIdx] = interiorDofIdx_;
}
else if (pressureDifference_[phaseIdx] < 0.0) {
upIdx_[phaseIdx] = interiorDofIdx_;
dnIdx_[phaseIdx] = exteriorDofIdx_;
}
else {
// if the pressure difference is zero, we chose the DOF which has the
// larger volume associated to it as upstream DOF
Scalar Vin = elemCtx.dofVolume(interiorDofIdx_, /*timeIdx=*/0);
Scalar Vex = elemCtx.dofVolume(exteriorDofIdx_, /*timeIdx=*/0);
if (Vin > Vex) {
upIdx_[phaseIdx] = interiorDofIdx_;
dnIdx_[phaseIdx] = exteriorDofIdx_;
}
else if (Vin < Vex) {
upIdx_[phaseIdx] = exteriorDofIdx_;
dnIdx_[phaseIdx] = interiorDofIdx_;
}
else {
assert(Vin == Vex);
// if the volumes are also equal, we pick the DOF which exhibits the
// smaller global index
if (I < J) {
upIdx_[phaseIdx] = interiorDofIdx_;
dnIdx_[phaseIdx] = exteriorDofIdx_;
}
else {
upIdx_[phaseIdx] = exteriorDofIdx_;
dnIdx_[phaseIdx] = interiorDofIdx_;
}
}
}
// apply the threshold pressure for the intersection. note that the concept
// of threshold pressure is a quite big hack that only makes sense for ECL
// datasets. (and even there, its physical justification is quite
// questionable IMO.)
if (std::abs(Toolbox::value(pressureDifference_[phaseIdx])) > thpres) {
if (pressureDifference_[phaseIdx] < 0.0)
pressureDifference_[phaseIdx] += thpres;
else
pressureDifference_[phaseIdx] -= thpres;
}
else {
pressureDifference_[phaseIdx] = 0.0;
volumeFlux_[phaseIdx] = 0.0;
continue;
}
// this is slightly hacky because in the automatic differentiation case, it
// only works for the element centered finite volume method. for ebos this
// does not matter, though.
unsigned upstreamIdx = upstreamIndex_(phaseIdx);
const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx);
if (upstreamIdx == interiorDofIdx_)
volumeFlux_[phaseIdx] =
pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-trans/faceArea);
else
volumeFlux_[phaseIdx] =
pressureDifference_[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*(-trans/faceArea));
}
}
/*!
* \brief Update the required gradients for boundary faces
*/
template <class FluidState>
void calculateBoundaryGradients_(const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx,
const FluidState& exFluidState)
{
const auto& problem = elemCtx.problem();
bool enableBoundaryMassFlux = problem.hasFreeBoundaryConditions();
if (!enableBoundaryMassFlux)
return;
const auto& stencil = elemCtx.stencil(timeIdx);
const auto& scvf = stencil.boundaryFace(scvfIdx);
interiorDofIdx_ = scvf.interiorIndex();
Scalar trans = problem.transmissibilityBoundary(elemCtx, scvfIdx);
Scalar faceArea = scvf.area();
// estimate the gravity correction: for performance reasons we use a simplified
// approach for this flux module that assumes that gravity is constant and always
// acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx_, timeIdx);
// this is quite hacky because the dune grid interface does not provide a
// cellCenterDepth() method (so we ask the problem to provide it). The "good"
// solution would be to take the Z coordinate of the element centroids, but since
// ECL seems to like to be inconsistent on that front, it needs to be done like
// here...
Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx_, timeIdx);
Scalar zEx = scvf.integrationPos()[dimWorld - 1];
// the distances from the DOF's depths. (i.e., the additional depth of the
// exterior DOF)
Scalar distZ = zIn - zEx;
for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
// do the gravity correction: compute the hydrostatic pressure for the
// integration position
const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
const auto& rhoEx = exFluidState.density(phaseIdx);
Evaluation rhoAvg = (rhoIn + rhoEx)/2;
const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
Evaluation pressureExterior = exFluidState.pressure(phaseIdx);
pressureExterior += rhoAvg*(distZ*g);
pressureDifference_[phaseIdx] = pressureExterior - pressureInterior;
// decide the upstream index for the phase. for this we make sure that the
// degree of freedom which is regarded upstream if both pressures are equal
// is always the same: if the pressure is equal, the DOF with the lower
// global index is regarded to be the upstream one.
if (pressureDifference_[phaseIdx] > 0.0) {
upIdx_[phaseIdx] = -1;
dnIdx_[phaseIdx] = interiorDofIdx_;
}
else {
upIdx_[phaseIdx] = interiorDofIdx_;
dnIdx_[phaseIdx] = -1;
}
// this is slightly hacky because in the automatic differentiation case, it
// only works for the element centered finite volume method. for ebos this
// does not matter, though.
unsigned upstreamIdx = upstreamIndex_(phaseIdx);
if (upstreamIdx == interiorDofIdx_) {
const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx);
volumeFlux_[phaseIdx] =
pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-trans/faceArea);
if (enableSolvent && phaseIdx == gasPhaseIdx) {
asImp_().setSolventVolumeFlux( pressureDifference_[phaseIdx]*up.solventMobility()*(-trans/faceArea));
}
}
else {
// compute the phase mobility using the material law parameters of the
// interior element. TODO: this could probably be done more efficiently
const auto& matParams =
elemCtx.problem().materialLawParams(elemCtx,
interiorDofIdx_,
/*timeIdx=*/0);
typename FluidState::Scalar kr[numPhases];
MaterialLaw::relativePermeabilities(kr, matParams, exFluidState);
const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx);
volumeFlux_[phaseIdx] =
pressureDifference_[phaseIdx]*mob*(-trans/faceArea);
// Solvent inflow is not yet supported
if (enableSolvent && phaseIdx == gasPhaseIdx) {
asImp_().setSolventVolumeFlux( 0.0 );
}
}
}
}
/*!
* \brief Update the volumetric fluxes for all fluid phases on the interior faces of the context
*/
void calculateFluxes_(const ElementContext& elemCtx OPM_UNUSED, unsigned scvfIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
{ }
void calculateBoundaryFluxes_(const ElementContext& elemCtx OPM_UNUSED, unsigned scvfIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
{}
private:
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
const Implementation& asImp_() const
{ return *static_cast<const Implementation*>(this); }
// the volumetric flux of all phases [m^3/s]
Evaluation volumeFlux_[numPhases];
// the difference in effective pressure between the exterior and the interior degree
// of freedom [Pa]
Evaluation pressureDifference_[numPhases];
// the local indices of the interior and exterior degrees of freedom
unsigned short interiorDofIdx_;
unsigned short exteriorDofIdx_;
unsigned short upIdx_[numPhases];
unsigned short dnIdx_[numPhases];
};
} // namespace Ewoms
#endif

232
ebos/eclnewtonmethod.hh Normal file
View File

@ -0,0 +1,232 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Ewoms::EclNewtonMethod
*/
#ifndef EWOMS_ECL_NEWTON_METHOD_HH
#define EWOMS_ECL_NEWTON_METHOD_HH
#include <ewoms/models/blackoil/blackoilnewtonmethod.hh>
#include <ewoms/common/signum.hh>
#include <opm/material/common/Unused.hpp>
BEGIN_PROPERTIES
NEW_PROP_TAG(EclNewtonSumTolerance);
NEW_PROP_TAG(EclNewtonStrictIterations);
NEW_PROP_TAG(EclNewtonRelaxedVolumeFraction);
NEW_PROP_TAG(EclNewtonSumToleranceExponent);
NEW_PROP_TAG(EclNewtonRelaxedTolerance);
END_PROPERTIES
namespace Ewoms {
/*!
* \brief A newton solver which is ebos specific.
*/
template <class TypeTag>
class EclNewtonMethod : public BlackOilNewtonMethod<TypeTag>
{
typedef BlackOilNewtonMethod<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) DiscNewtonMethod;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Linearizer) Linearizer;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
static const unsigned numEq = GET_PROP_VALUE(TypeTag, NumEq);
static constexpr int contiSolventEqIdx = Indices::contiSolventEqIdx;
static constexpr int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
static constexpr int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
friend NewtonMethod<TypeTag>;
friend DiscNewtonMethod;
friend ParentType;
public:
EclNewtonMethod(Simulator& simulator) : ParentType(simulator)
{
errorPvFraction_ = 1.0;
relaxedMaxPvFraction_ = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonRelaxedVolumeFraction);
sumTolerance_ = 0.0; // this gets determined in the error calculation proceedure
relaxedTolerance_ = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonRelaxedTolerance);
numStrictIterations_ = EWOMS_GET_PARAM(TypeTag, int, EclNewtonStrictIterations);
}
/*!
* \brief Register all run-time parameters for the Newton method.
*/
static void registerParameters()
{
ParentType::registerParameters();
EWOMS_REGISTER_PARAM(TypeTag, Scalar, EclNewtonSumTolerance,
"The maximum error tolerated by the Newton"
"method for considering a solution to be "
"converged");
EWOMS_REGISTER_PARAM(TypeTag, int, EclNewtonStrictIterations,
"The number of Newton iterations where the"
" volumetric error is considered.");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, EclNewtonRelaxedVolumeFraction,
"The fraction of the pore volume of the reservoir "
"where the volumetric error may be voilated during "
"strict Newton iterations.");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, EclNewtonSumToleranceExponent,
"The the exponent used to scale the sum tolerance by "
"the total pore volume of the reservoir.");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, EclNewtonRelaxedTolerance,
"The maximum error which the volumetric residual "
"may exhibit if it is in a 'relaxed' "
"region during a strict iteration.");
}
/*!
* \brief Returns true if the error of the solution is below the
* tolerance.
*/
bool converged() const
{
if (errorPvFraction_ < relaxedMaxPvFraction_)
return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ;
else if (this->numIterations() > numStrictIterations_)
return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ;
return this->error_ <= this->tolerance() && errorSum_ <= sumTolerance_;
}
void preSolve_(const SolutionVector& currentSolution OPM_UNUSED,
const GlobalEqVector& currentResidual)
{
const auto& constraintsMap = this->model().linearizer().constraintsMap();
this->lastError_ = this->error_;
Scalar newtonMaxError = EWOMS_GET_PARAM(TypeTag, Scalar, NewtonMaxError);
// calculate the error as the maximum weighted tolerance of
// the solution's residual
this->error_ = 0.0;
Dune::FieldVector<Scalar, numEq> componentSumError;
std::fill(componentSumError.begin(), componentSumError.end(), 0.0);
Scalar sumPv = 0.0;
errorPvFraction_ = 0.0;
const Scalar dt = this->simulator_.timeStepSize();
for (unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) {
// do not consider auxiliary DOFs for the error
if (dofIdx >= this->model().numGridDof()
|| this->model().dofTotalVolume(dofIdx) <= 0.0)
continue;
if (!this->model().isLocalDof(dofIdx))
continue;
// also do not consider DOFs which are constraint
if (this->enableConstraints_()) {
if (constraintsMap.count(dofIdx) > 0)
continue;
}
const auto& r = currentResidual[dofIdx];
const double pvValue =
this->simulator_.problem().porosity(dofIdx)
* this->model().dofTotalVolume(dofIdx);
sumPv += pvValue;
bool cnvViolated = false;
for (unsigned eqIdx = 0; eqIdx < r.size(); ++eqIdx) {
Scalar tmpError = r[eqIdx] * dt * this->model().eqWeight(dofIdx, eqIdx) / pvValue;
Scalar tmpError2 = r[eqIdx] * this->model().eqWeight(dofIdx, eqIdx);
this->error_ = Opm::max(std::abs(tmpError), this->error_);
if (std::abs(tmpError) > this->tolerance_)
cnvViolated = true;
componentSumError[eqIdx] += std::abs(tmpError2);
}
if (cnvViolated)
errorPvFraction_ += pvValue;
}
// take the other processes into account
this->error_ = this->comm_.max(this->error_);
componentSumError = this->comm_.sum(componentSumError);
sumPv = this->comm_.sum(sumPv);
errorPvFraction_ = this->comm_.sum(errorPvFraction_);
componentSumError /= sumPv;
componentSumError *= dt;
errorPvFraction_ /= sumPv;
errorSum_ = 0;
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
errorSum_ = std::max(std::abs(componentSumError[eqIdx]), errorSum_);
// scale the tolerance for the total error with the pore volume. by default, the
// exponent is 1/3, i.e., cubic root.
Scalar x = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonSumTolerance);
Scalar y = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonSumToleranceExponent);
sumTolerance_ = x*std::pow(sumPv, y);
// make sure that the error never grows beyond the maximum
// allowed one
if (this->error_ > newtonMaxError)
throw Opm::NumericalIssue("Newton: Error "+std::to_string(double(this->error_))
+" is larger than maximum allowed error of "
+std::to_string(double(newtonMaxError)));
// make sure that the error never grows beyond the maximum
// allowed one
if (errorSum_ > newtonMaxError)
throw Opm::NumericalIssue("Newton: Sum of the error "+std::to_string(double(errorSum_))
+" is larger than maximum allowed error of "
+std::to_string(double(newtonMaxError)));
}
private:
Scalar errorPvFraction_;
Scalar errorSum_;
Scalar relaxedTolerance_;
Scalar relaxedMaxPvFraction_;
Scalar sumTolerance_;
int numStrictIterations_;
};
} // namespace Ewoms
#endif

File diff suppressed because it is too large Load Diff

1722
ebos/eclpeacemanwell.hh Normal file

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,170 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
* \copydoc Ewoms::EclPolyhedralGridVanguard
*/
#ifndef EWOMS_ECL_POLYHEDRAL_GRID_VANGUARD_HH
#define EWOMS_ECL_POLYHEDRAL_GRID_VANGUARD_HH
#include "eclbasevanguard.hh"
#include <opm/grid/polyhedralgrid.hh>
namespace Ewoms {
template <class TypeTag>
class EclPolyhedralGridVanguard;
}
BEGIN_PROPERTIES
NEW_TYPE_TAG(EclPolyhedralGridVanguard, INHERITS_FROM(EclBaseVanguard));
// declare the properties
SET_TYPE_PROP(EclPolyhedralGridVanguard, Vanguard, Ewoms::EclPolyhedralGridVanguard<TypeTag>);
SET_TYPE_PROP(EclPolyhedralGridVanguard, Grid, Dune::PolyhedralGrid<3, 3>);
SET_TYPE_PROP(EclPolyhedralGridVanguard, EquilGrid, typename GET_PROP_TYPE(TypeTag, Grid));
END_PROPERTIES
namespace Ewoms {
/*!
* \ingroup EclBlackOilSimulator
*
* \brief Helper class for grid instantiation of ECL file-format using problems.
*
* This class uses Dune::PolyhedralGrid as the simulation grid.
*/
template <class TypeTag>
class EclPolyhedralGridVanguard : public EclBaseVanguard<TypeTag>
{
friend class EclBaseVanguard<TypeTag>;
typedef EclBaseVanguard<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
public:
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, EquilGrid) EquilGrid;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
private:
typedef Grid* GridPointer;
typedef EquilGrid* EquilGridPointer;
typedef Dune::CartesianIndexMapper<Grid> CartesianIndexMapper;
typedef CartesianIndexMapper* CartesianIndexMapperPointer;
public:
/*!
* \brief Inherit the constructors from the base class.
*/
using EclBaseVanguard<TypeTag>::EclBaseVanguard;
~EclPolyhedralGridVanguard()
{
delete cartesianIndexMapper_;
delete grid_;
}
/*!
* \brief Return a reference to the simulation grid.
*/
Grid& grid()
{ return *grid_; }
/*!
* \brief Return a reference to the simulation grid.
*/
const Grid& grid() const
{ return *grid_; }
/*!
* \brief Returns a refefence to the grid which should be used by the EQUIL
* initialization code.
*
* The EQUIL keyword is used to specify the initial condition of the reservoir in
* hydrostatic equilibrium. Since the code which does this is not accepting arbitrary
* DUNE grids (the code is part of the opm-core module), this is not necessarily the
* same as the grid which is used for the actual simulation.
*/
const EquilGrid& equilGrid() const
{ return *grid_; }
/*!
* \brief Indicates that the initial condition has been computed and the memory used
* by the EQUIL grid can be released.
*
* Depending on the implementation, subsequent accesses to the EQUIL grid lead to
* crashes.
*/
void releaseEquilGrid()
{ /* do nothing: The EQUIL grid is the simulation grid! */ }
/*!
* \brief Distribute the simulation grid over multiple processes
*
* (For parallel simulation runs.)
*/
void loadBalance()
{ /* do nothing: PolyhedralGrid is not parallel! */ }
/*!
* \brief Returns the object which maps a global element index of the simulation grid
* to the corresponding element index of the logically Cartesian index.
*/
const CartesianIndexMapper& cartesianIndexMapper() const
{ return *cartesianIndexMapper_; }
/*!
* \brief Returns mapper from compressed to cartesian indices for the EQUIL grid
*
* Since PolyhedralGrid is not parallel, that's always the same as
* cartesianIndexMapper().
*/
const CartesianIndexMapper& equilCartesianIndexMapper() const
{ return *cartesianIndexMapper_; }
protected:
void createGrids_()
{
const auto& gridProps = this->eclState().get3DProperties();
const std::vector<double>& porv = gridProps.getDoubleGridProperty("PORV").getData();
grid_ = new Grid(this->deck(), porv);
cartesianIndexMapper_ = new CartesianIndexMapper(*grid_);
}
void filterConnections_()
{
// not handling the removal of completions for this type of grid yet.
}
GridPointer grid_;
CartesianIndexMapperPointer cartesianIndexMapper_;
};
} // namespace Ewoms
#endif

2642
ebos/eclproblem.hh Normal file

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,324 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Ewoms::EclThresholdPressure
*/
#ifndef EWOMS_ECL_THRESHOLD_PRESSURE_HH
#define EWOMS_ECL_THRESHOLD_PRESSURE_HH
#include <ewoms/common/propertysystem.hh>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/Eqldims.hpp>
#include <opm/parser/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
#include <opm/parser/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <dune/grid/common/gridenums.hh>
#include <dune/common/version.hh>
#include <array>
#include <vector>
#include <unordered_map>
BEGIN_PROPERTIES
NEW_PROP_TAG(Simulator);
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(Evaluation);
NEW_PROP_TAG(ElementContext);
NEW_PROP_TAG(FluidSystem);
END_PROPERTIES
namespace Ewoms {
/*!
* \ingroup EclBlackOilSimulator
*
* \brief This class calculates the threshold pressure for grid faces according to the
* Eclipse Reference Manual.
*
* If the difference of the pressure potential between two cells is below the threshold
* pressure, the pressure potential difference is assumed to be zero, if it is larger
* than the threshold pressure, it is reduced by the threshold pressure.
*/
template <class TypeTag>
class EclThresholdPressure
{
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
enum { numPhases = FluidSystem::numPhases };
public:
EclThresholdPressure(const Simulator& simulator)
: simulator_(simulator)
{
enableThresholdPressure_ = false;
}
void setFromRestart(const std::vector<Scalar>& values)
{ thpres_ = values; }
/*!
* \brief Actually compute the threshold pressures over a face as a pre-compute step.
*/
void finishInit()
{
const auto& gridView = simulator_.gridView();
unsigned numElements = gridView.size(/*codim=*/0);
// this code assumes that the DOFs are the elements. (i.e., an
// ECFV spatial discretization with TPFA). if you try to use
// it with something else, you're currently out of luck,
// sorry!
assert(simulator_.model().numGridDof() == numElements);
const auto& vanguard = simulator_.vanguard();
const auto& eclState = vanguard.eclState();
const auto& simConfig = eclState.getSimulationConfig();
enableThresholdPressure_ = simConfig.useThresholdPressure();
if (!enableThresholdPressure_)
return;
numEquilRegions_ = eclState.getTableManager().getEqldims().getNumEquilRegions();
if (numEquilRegions_ > 0xff) {
// make sure that the index of an equilibration region can be stored in a
// single byte
throw std::runtime_error("The maximum number of supported equilibration regions is 255!");
}
// internalize the data specified using the EQLNUM keyword
const std::vector<int>& equilRegionData =
eclState.get3DProperties().getIntGridProperty("EQLNUM").getData();
elemEquilRegion_.resize(numElements, 0);
for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
int cartElemIdx = vanguard.cartesianIndex(elemIdx);
// ECL uses Fortran-style indices but we want C-style ones!
elemEquilRegion_[elemIdx] = equilRegionData[cartElemIdx] - 1;
}
/*
If this is a restart run the ThresholdPressure object will be active,
but it will *not* be properly initialized with numerical values. The
values must instead come from the THPRES vector in the restart file.
*/
if (simConfig.getThresholdPressure().restart())
return;
// allocate the array which specifies the threshold pressures
thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0);
thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0);
computeDefaultThresholdPressures_();
applyExplicitThresholdPressures_();
}
/*!
* \brief Returns the theshold pressure [Pa] for the intersection between two elements.
*
* This is tailor made for the E100 threshold pressure mechanism and it is thus quite
* a hack: First of all threshold pressures in general are unphysical, and second,
* they should be different for the fluid phase but are not. Anyway, this seems to be
* E100's way of doing things, so we do it the same way.
*/
Scalar thresholdPressure(int elemIdx1, int elemIdx2) const
{
if (!enableThresholdPressure_)
return 0.0;
unsigned short equilRegion1Idx = elemEquilRegion_[elemIdx1];
unsigned short equilRegion2Idx = elemEquilRegion_[elemIdx2];
if (equilRegion1Idx == equilRegion2Idx)
return 0.0;
return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx];
}
const std::vector<Scalar>& data() const {
return thpres_;
}
private:
// compute the defaults of the threshold pressures using the initial condition
void computeDefaultThresholdPressures_()
{
const auto& vanguard = simulator_.vanguard();
const auto& gridView = vanguard.gridView();
typedef Opm::MathToolbox<Evaluation> Toolbox;
// loop over the whole grid and compute the maximum gravity adjusted pressure
// difference between two EQUIL regions.
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
ElementContext elemCtx(simulator_);
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity)
continue;
elemCtx.updateAll(elem);
const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
for (unsigned scvfIdx = 0; scvfIdx < stencil.numInteriorFaces(); ++ scvfIdx) {
const auto& face = stencil.interiorFace(scvfIdx);
unsigned i = face.interiorIndex();
unsigned j = face.exteriorIndex();
unsigned insideElemIdx = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0);
unsigned outsideElemIdx = elemCtx.globalSpaceIndex(j, /*timeIdx=*/0);
unsigned equilRegionInside = elemEquilRegion_[insideElemIdx];
unsigned equilRegionOutside = elemEquilRegion_[outsideElemIdx];
if (equilRegionInside == equilRegionOutside)
// the current face is not at the boundary between EQUIL regions!
continue;
// don't include connections with negligible flow
const Scalar& trans = simulator_.problem().transmissibility(elemCtx, i, j);
const Scalar& faceArea = face.area();
if ( std::abs(faceArea * trans) < 1e-18)
continue;
// determine the maximum difference of the pressure of any phase over the
// intersection
Scalar pth = 0.0;
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx=*/0);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
const auto& up = elemCtx.intensiveQuantities(upIdx, /*timeIdx=*/0);
if (up.mobility(phaseIdx) > 0.0) {
Scalar phaseVal = Toolbox::value(extQuants.pressureDifference(phaseIdx));
pth = std::max(pth, std::abs(phaseVal));
}
}
int offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside;
int offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside;
thpresDefault_[offset1] = std::max(thpresDefault_[offset1], pth);
thpresDefault_[offset2] = std::max(thpresDefault_[offset2], pth);
}
}
// make sure that the threshold pressures is consistent for parallel
// runs. (i.e. take the maximum of all processes)
for (unsigned i = 0; i < thpresDefault_.size(); ++i)
thpresDefault_[i] = gridView.comm().max(thpresDefault_[i]);
}
// internalize the threshold pressures which where explicitly specified via the
// THPRES keyword.
void applyExplicitThresholdPressures_()
{
const auto& vanguard = simulator_.vanguard();
const auto& gridView = vanguard.gridView();
const auto& elementMapper = simulator_.model().elementMapper();
const auto& eclState = simulator_.vanguard().eclState();
const Opm::SimulationConfig& simConfig = eclState.getSimulationConfig();
const auto& thpres = simConfig.getThresholdPressure();
// set the threshold pressures for all EQUIL region boundaries which have a
// intersection in the grid
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity)
continue;
auto isIt = gridView.ibegin(elem);
const auto& isEndIt = gridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
// store intersection, this might be costly
const auto& intersection = *isIt;
// ignore boundary intersections for now (TODO?)
if (intersection.boundary())
continue;
const auto& inside = intersection.inside();
const auto& outside = intersection.outside();
unsigned insideElemIdx = elementMapper.index(inside);
unsigned outsideElemIdx = elementMapper.index(outside);
unsigned equilRegionInside = elemEquilRegion_[insideElemIdx];
unsigned equilRegionOutside = elemEquilRegion_[outsideElemIdx];
if (thpres.hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1)) {
Scalar pth = 0.0;
if (thpres.hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) {
// threshold pressure explicitly specified
pth = thpres.getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1);
}
else {
// take the threshold pressure from the initial condition
unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside;
pth = thpresDefault_[offset];
}
unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside;
unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside;
thpres_[offset1] = pth;
thpres_[offset2] = pth;
}
}
}
}
const Simulator& simulator_;
std::vector<Scalar> thpresDefault_;
std::vector<Scalar> thpres_;
unsigned numEquilRegions_;
std::vector<unsigned char> elemEquilRegion_;
bool enableThresholdPressure_;
};
} // namespace Ewoms
#endif

517
ebos/ecltracermodel.hh Normal file
View File

@ -0,0 +1,517 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/**
* \file
*
* \copydoc Ewoms::EclTracerModel
*/
#ifndef EWOMS_ECL_TRACER_MODEL_HH
#define EWOMS_ECL_TRACER_MODEL_HH
#include "tracervdtable.hh"
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
#include <string>
#include <vector>
#include <iostream>
BEGIN_PROPERTIES
NEW_PROP_TAG(EnableTracerModel);
END_PROPERTIES
namespace Ewoms {
/*!
* \ingroup EclBlackOilSimulator
*
* \brief A class which handles tracers as specified in by ECL
*
* TODO: MPI parallelism.
*/
template <class TypeTag>
class EclTracerModel
{
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Stencil) Stencil;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef Opm::DenseAd::Evaluation<Scalar,1> TracerEvaluation;
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
enum { numPhases = FluidSystem::numPhases };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 1, 1>> TracerMatrix;
typedef Dune::BlockVector<Dune::FieldVector<Scalar,1>> TracerVector;
public:
EclTracerModel(Simulator& simulator)
: simulator_(simulator)
{ }
/*!
* \brief Initialize all internal data structures needed by the tracer module
*/
void init()
{
const Opm::Deck& deck = simulator_.vanguard().deck();
if (!deck.hasKeyword("TRACERS"))
return; // tracer treatment is supposed to be disabled
if (!EWOMS_GET_PARAM(TypeTag, bool, EnableTracerModel)) {
std::cout << "Warning: Tracer model is disabled but the deck contatins the TRACERS keyword \n";
std::cout << "The tracer model must be activated using --enable-tracer-model=true "<< std::endl;
return; // Tracer transport must be enabled by the user
}
if (!deck.hasKeyword("TRACER"))
throw std::runtime_error("the deck does not contain the TRACER keyword");
if (simulator_.gridView().comm().size() > 1) {
tracerNames_.resize(0);
std::cout << "Warning: Tracer model is not compatible with mpi run" << std::endl;
return;
}
// retrieve the number of tracers from the deck
const int numTracers = deck.getKeyword("TRACER").size();
tracerNames_.resize(numTracers);
tracerConcentration_.resize(numTracers);
storageOfTimeIndex1_.resize(numTracers);
// the phase where the tracer is
tracerPhaseIdx_.resize(numTracers);
size_t numGridDof = simulator_.model().numGridDof();
for (int tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx) {
const auto& tracerRecord = deck.getKeyword("TRACER").getRecord(tracerIdx);
tracerNames_[tracerIdx] = tracerRecord.getItem("NAME").template get<std::string>(0);
const std::string& fluidName = tracerRecord.getItem("FLUID").template get<std::string>(0);
if (fluidName == "WAT")
tracerPhaseIdx_[tracerIdx] = waterPhaseIdx;
else if (fluidName == "OIL")
tracerPhaseIdx_[tracerIdx] = oilPhaseIdx;
else if (fluidName == "GAS")
tracerPhaseIdx_[tracerIdx] = gasPhaseIdx;
else
throw std::invalid_argument("Tracer: invalid fluid name "
+fluidName+" for "+tracerNames_[tracerIdx]);
tracerConcentration_[tracerIdx].resize(numGridDof);
storageOfTimeIndex1_[tracerIdx].resize(numGridDof);
std::string tmp = "TVDPF" +tracerNames_[tracerIdx];
//TBLK keyword
if (deck.hasKeyword("TBLKF" +tracerNames_[tracerIdx])){
const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper();
const auto& tblkData =
deck.getKeyword("TBLKF"
+tracerNames_
[tracerIdx]).getRecord(0).getItem(0).getSIDoubleData();
int tblkDatasize = tblkData.size();
if (tblkDatasize < simulator_.vanguard().cartesianSize()){
throw std::runtime_error("Uninitialized tracer concentration (TBLKF) for tracer "
+ tracerName(tracerIdx));
}
for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
int cartDofIdx = cartMapper.cartesianIndex(globalDofIdx);
tracerConcentration_[tracerIdx][globalDofIdx] = tblkData[cartDofIdx];
}
}
//TVDPF keyword
else if (deck.hasKeyword(tmp)){
TracerVdTable dtable(deck.getKeyword(tmp).getRecord(0).getItem(0));
const auto& eclGrid = simulator_.vanguard().eclState().getInputGrid();
const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper();
for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
int cartDofIdx = cartMapper.cartesianIndex(globalDofIdx);
const auto& center = eclGrid.getCellCenter(cartDofIdx);
tracerConcentration_[tracerIdx][globalDofIdx]
= dtable.evaluate("TRACER_CONCENTRATION", center[2]);
}
}
else {
throw std::runtime_error("Uninitialized tracer concentration for tracer "
+ tracerName(tracerIdx));
}
}
// initial tracer concentration
tracerConcentrationInitial_ = tracerConcentration_;
// residual of tracers
tracerResidual_.resize(numGridDof);
// allocate matrix for storing the Jacobian of the tracer residual
tracerMatrix_ = new TracerMatrix(numGridDof, numGridDof, TracerMatrix::random);
// find the sparsity pattern of the tracer matrix
typedef std::set<unsigned> NeighborSet;
std::vector<NeighborSet> neighbors(numGridDof);
Stencil stencil(simulator_.gridView(), simulator_.model().dofMapper() );
ElementIterator elemIt = simulator_.gridView().template begin<0>();
const ElementIterator elemEndIt = simulator_.gridView().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
stencil.update(elem);
for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
neighbors[myIdx].insert(neighborIdx);
}
}
}
// allocate space for the rows of the matrix
for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx)
tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
tracerMatrix_->endrowsizes();
// fill the rows with indices. each degree of freedom talks to
// all of its neighbors. (it also talks to itself since
// degrees of freedom are sometimes quite egocentric.)
for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
for (; nIt != nEndIt; ++nIt)
tracerMatrix_->addindex(dofIdx, *nIt);
}
tracerMatrix_->endindices();
const int sizeCartGrid = simulator_.vanguard().cartesianSize();
cartToGlobal_.resize(sizeCartGrid);
for (unsigned i = 0; i < numGridDof; ++i) {
int cartIdx = simulator_.vanguard().cartesianIndex(i);
cartToGlobal_[cartIdx] = i;
}
}
/*!
* \brief Return the number of tracers considered by the tracerModel.
*/
int numTracers() const
{ return tracerNames_.size(); }
/*!
* \brief Return the tracer name
*/
const std::string& tracerName(int tracerIdx) const
{
if (numTracers()==0)
throw std::logic_error("This method should never be called when there are no tracers in the model");
return tracerNames_[tracerIdx];
}
/*!
* \brief Return the tracer concentration for tracer index and global DofIdx
*/
Scalar tracerConcentration(int tracerIdx, int globalDofIdx) const
{
if (numTracers()==0)
return 0.0;
return tracerConcentration_[tracerIdx][globalDofIdx];
}
void beginTimeStep()
{
if (numTracers()==0)
return;
tracerConcentrationInitial_ = tracerConcentration_;
// compute storageCache
ElementContext elemCtx(simulator_);
auto elemIt = simulator_.gridView().template begin</*codim=*/0>();
auto elemEndIt = simulator_.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
elemCtx.updateAll(*elemIt);
int globalDofIdx = elemCtx.globalSpaceIndex(0, 0);
for (int tracerIdx = 0; tracerIdx < numTracers(); ++ tracerIdx){
Scalar storageOfTimeIndex1;
computeStorage_(storageOfTimeIndex1, elemCtx, 0, /*timIdx=*/0, tracerIdx);
storageOfTimeIndex1_[tracerIdx][globalDofIdx] = storageOfTimeIndex1;
}
}
}
/*!
* \brief Informs the tracer model that a time step has just been finished.
*/
void endTimeStep()
{
if (numTracers()==0)
return;
for (int tracerIdx = 0; tracerIdx < numTracers(); ++ tracerIdx){
TracerVector dx(tracerResidual_.size());
// Newton step (currently the system is linear, converge in one iteration)
for (int iter = 0; iter < 5; ++ iter){
linearize_(tracerIdx);
linearSolve_(*tracerMatrix_, dx, tracerResidual_);
tracerConcentration_[tracerIdx] -= dx;
if (dx.two_norm()<1e-2)
break;
}
}
}
/*!
* \brief This method writes the complete state of all tracer
* to the hard disk.
*/
template <class Restarter>
void serialize(Restarter& res OPM_UNUSED)
{ /* not implemented */ }
/*!
* \brief This method restores the complete state of the tracer
* from disk.
*
* It is the inverse of the serialize() method.
*/
template <class Restarter>
void deserialize(Restarter& res OPM_UNUSED)
{ /* not implemented */ }
protected:
// evaluate storage term for all tracers in a single cell
template <class LhsEval>
void computeStorage_(LhsEval& tracerStorage,
const ElementContext& elemCtx,
unsigned scvIdx,
unsigned timeIdx,
const int tracerIdx)
{
int globalDofIdx = elemCtx.globalSpaceIndex(scvIdx, timeIdx);
const auto& intQuants = elemCtx.intensiveQuantities(scvIdx, timeIdx);
const auto& fs = intQuants.fluidState();
Scalar phaseVolume =
Opm::decay<Scalar>(fs.saturation(tracerPhaseIdx_[tracerIdx]))
*Opm::decay<Scalar>(fs.invB(tracerPhaseIdx_[tracerIdx]))
*Opm::decay<Scalar>(intQuants.porosity());
// avoid singular matrix if no water is present.
phaseVolume = Opm::max(phaseVolume, 1e-10);
if (std::is_same<LhsEval, Scalar>::value)
tracerStorage = phaseVolume * tracerConcentrationInitial_[tracerIdx][globalDofIdx];
else
tracerStorage =
phaseVolume
* Opm::variable<LhsEval>(tracerConcentration_[tracerIdx][globalDofIdx][0], 0);
}
// evaluate the tracerflux over one face
void computeFlux_(TracerEvaluation & tracerFlux,
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx,
const int tracerIdx)
{
const auto& stencil = elemCtx.stencil(timeIdx);
const auto& scvf = stencil.interiorFace(scvfIdx);
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
unsigned inIdx = extQuants.interiorIndex();
const int tracerPhaseIdx = tracerPhaseIdx_[tracerIdx];
unsigned upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
int globalUpIdx = elemCtx.globalSpaceIndex(upIdx, timeIdx);
const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
const auto& fs = intQuants.fluidState();
Scalar A = scvf.area();
Scalar v = Opm::decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx));
Scalar b = Opm::decay<Scalar>(fs.invB(tracerPhaseIdx_[tracerIdx]));
Scalar c = tracerConcentration_[tracerIdx][globalUpIdx];
if (inIdx == upIdx)
tracerFlux = A*v*b*Opm::variable<TracerEvaluation>(c, 0);
else
tracerFlux = A*v*b*c;
}
bool linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
{
#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
Dune::FMatrixPrecision<Scalar>::set_singular_limit(1.e-30);
Dune::FMatrixPrecision<Scalar>::set_absolute_limit(1.e-30);
#endif
x = 0.0;
Scalar tolerance = 1e-2;
int maxIter = 100;
int verbosity = 0;
typedef Dune::BiCGSTABSolver<TracerVector> TracerSolver;
typedef Dune::MatrixAdapter<TracerMatrix, TracerVector , TracerVector > TracerOperator;
typedef Dune::SeqScalarProduct< TracerVector > TracerScalarProduct ;
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2,6)
typedef Dune::SeqILU< TracerMatrix, TracerVector, TracerVector > TracerPreconditioner;
#else
typedef Dune::SeqILUn< TracerMatrix, TracerVector, TracerVector > TracerPreconditioner;
#endif
TracerOperator tracerOperator(M);
TracerScalarProduct tracerScalarProduct;
TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
TracerSolver solver (tracerOperator, tracerScalarProduct,
tracerPreconditioner, tolerance, maxIter,
verbosity);
Dune::InverseOperatorResult result;
solver.apply(x, b, result);
// return the result of the solver
return result.converged;
}
void linearize_(int tracerIdx)
{
(*tracerMatrix_) = 0.0;
tracerResidual_ = 0.0;
size_t numGridDof = simulator_.model().numGridDof();
std::vector<double> volumes(numGridDof, 0.0);
ElementContext elemCtx(simulator_);
auto elemIt = simulator_.gridView().template begin</*codim=*/0>();
auto elemEndIt = simulator_.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
elemCtx.updateAll(*elemIt);
Scalar extrusionFactor =
elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
Opm::Valgrind::CheckDefined(extrusionFactor);
assert(Opm::isfinite(extrusionFactor));
assert(extrusionFactor > 0.0);
Scalar scvVolume =
elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume()
* extrusionFactor;
Scalar dt = elemCtx.simulator().timeStepSize();
size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timIdx=*/0);
volumes[I] = scvVolume;
TracerEvaluation localStorage;
TracerEvaluation storageOfTimeIndex0;
Scalar storageOfTimeIndex1;
computeStorage_(storageOfTimeIndex0, elemCtx, 0, /*timIdx=*/0, tracerIdx);
if (elemCtx.enableStorageCache())
storageOfTimeIndex1 = storageOfTimeIndex1_[tracerIdx][I];
else
computeStorage_(storageOfTimeIndex1, elemCtx, 0, /*timIdx=*/1, tracerIdx);
localStorage = (storageOfTimeIndex0 - storageOfTimeIndex1) * scvVolume/dt;
tracerResidual_[I][0] += localStorage.value(); //residual + flux
(*tracerMatrix_)[I][I][0][0] = localStorage.derivative(0);
size_t numInteriorFaces = elemCtx.numInteriorFaces(/*timIdx=*/0);
for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
TracerEvaluation flux;
const auto& face = elemCtx.stencil(0).interiorFace(scvfIdx);
unsigned j = face.exteriorIndex();
unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
computeFlux_(flux, elemCtx, scvfIdx, 0, tracerIdx);
tracerResidual_[I][0] += flux.value(); //residual + flux
(*tracerMatrix_)[J][I][0][0] = -flux.derivative(0);
(*tracerMatrix_)[I][J][0][0] = flux.derivative(0);
}
}
// Wells
const int episodeIdx = simulator_.episodeIndex();
const auto& wells = simulator_.vanguard().schedule().getWells(episodeIdx);
for (auto well : wells) {
if (well->getStatus(episodeIdx) == Opm::WellCommon::SHUT)
continue;
const double wtracer = well->getTracerProperties(episodeIdx).getConcentration(tracerNames_[tracerIdx]);
std::array<int, 3> cartesianCoordinate;
for (auto& connection : well->getConnections(episodeIdx)) {
if (connection.state() == Opm::WellCompletion::SHUT)
continue;
cartesianCoordinate[0] = connection.getI();
cartesianCoordinate[1] = connection.getJ();
cartesianCoordinate[2] = connection.getK();
const size_t cartIdx = simulator_.vanguard().cartesianIndex(cartesianCoordinate);
const int I = cartToGlobal_[cartIdx];
Scalar rate = simulator_.problem().wellModel().well(well->name())->volumetricSurfaceRateForConnection(I, tracerPhaseIdx_[tracerIdx]);
if (rate > 0)
tracerResidual_[I][0] -= rate*wtracer;
else if (rate < 0)
tracerResidual_[I][0] -= rate*tracerConcentration_[tracerIdx][I];
}
}
}
Simulator& simulator_;
std::vector<std::string> tracerNames_;
std::vector<int> tracerPhaseIdx_;
std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> tracerConcentration_;
std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> tracerConcentrationInitial_;
TracerMatrix *tracerMatrix_;
TracerVector tracerResidual_;
std::vector<int> cartToGlobal_;
std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> storageOfTimeIndex1_;
};
} // namespace Ewoms
#endif

953
ebos/ecltransmissibility.hh Normal file
View File

@ -0,0 +1,953 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Ewoms::EclTransmissibility
*/
#ifndef EWOMS_ECL_TRANSMISSIBILITY_HH
#define EWOMS_ECL_TRANSMISSIBILITY_HH
#include <ewoms/common/propertysystem.hh>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/GridProperties.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/TransMult.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <opm/material/common/ConditionalStorage.hpp>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/common/version.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <array>
#include <vector>
#include <unordered_map>
BEGIN_PROPERTIES
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(Vanguard);
NEW_PROP_TAG(Grid);
NEW_PROP_TAG(GridView);
NEW_PROP_TAG(ElementMapper);
NEW_PROP_TAG(EnableEnergy);
END_PROPERTIES
namespace Ewoms {
/*!
* \ingroup EclBlackOilSimulator
*
* \brief This class calculates the transmissibilites for grid faces according to the
* Eclipse Technical Description.
*/
template <class TypeTag>
class EclTransmissibility
{
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard;
typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
typedef typename GridView::Intersection Intersection;
static const bool enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy);
// Grid and world dimension
enum { dimWorld = GridView::dimensionworld };
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
static const unsigned elemIdxShift = 32; // bits
public:
EclTransmissibility(const Vanguard& vanguard)
: vanguard_(vanguard)
{
const Opm::UnitSystem& unitSystem = vanguard_.deck().getActiveUnitSystem();
transmissibilityThreshold_ = unitSystem.parse("Transmissibility").getSIScaling() * 1e-6;
}
/*!
* \brief Actually compute the transmissibilty over a face as a pre-compute step.
*
* This code actually uses the direction specific "centroids" of
* each element. These "centroids" are _not_ the identical
* barycenter of the element, but the middle of the centers of the
* faces of the logical Cartesian cells, i.e., the centers of the
* faces of the reference elements. We do things this way because
* the barycenter of the element can be located outside of the
* element for sufficiently "ugly" (i.e., thin and "non-flat")
* elements which in turn leads to quite wrong
* permeabilities. This approach is probably not always correct
* either but at least it seems to be much better.
*/
void finishInit()
{ update(); }
/*!
* \brief Compute all transmissibilities
*
* Also, this updates the "thermal half transmissibilities" if energy is enabled.
*/
void update()
{
const auto& gridView = vanguard_.gridView();
const auto& cartMapper = vanguard_.cartesianIndexMapper();
const auto& eclState = vanguard_.eclState();
const auto& eclGrid = eclState.getInputGrid();
const auto& cartDims = cartMapper.cartesianDimensions();
auto& transMult = eclState.getTransMult();
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
#else
ElementMapper elemMapper(gridView);
#endif
// get the ntg values, the ntg values are modified for the cells merged with minpv
std::vector<double> ntg;
minPvFillNtg_(ntg);
unsigned numElements = elemMapper.size();
extractPermeability_();
// calculate the axis specific centroids of all elements
std::array<std::vector<DimVector>, dimWorld> axisCentroids;
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
axisCentroids[dimIdx].resize(numElements);
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
unsigned elemIdx = elemMapper.index(elem);
// compute the axis specific "centroids" used for the transmissibilities. for
// consistency with the flow simulator, we use the element centers as
// computed by opm-parser's Opm::EclipseGrid class for all axes.
unsigned cartesianCellIdx = cartMapper.cartesianIndex(elemIdx);
const auto& centroid = eclGrid.getCellCenter(cartesianCellIdx);
for (unsigned axisIdx = 0; axisIdx < dimWorld; ++axisIdx)
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
axisCentroids[axisIdx][elemIdx][dimIdx] = centroid[dimIdx];
}
// reserving some space in the hashmap upfront saves quite a bit of time because
// resizes are costly for hashmaps and there would be quite a few of them if we
// would not have a rough idea of how large the final map will be (the rough idea
// is a conforming Cartesian grid).
trans_.clear();
trans_.reserve(numElements*3*1.05);
transBoundary_.clear();
// if energy is enabled, let's do the same for the "thermal half transmissibilities"
if (enableEnergy) {
thermalHalfTrans_->clear();
thermalHalfTrans_->reserve(numElements*6*1.05);
thermalHalfTransBoundary_.clear();
}
// compute the transmissibilities for all intersections
elemIt = gridView.template begin</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
unsigned elemIdx = elemMapper.index(elem);
auto isIt = gridView.ibegin(elem);
const auto& isEndIt = gridView.iend(elem);
unsigned boundaryIsIdx = 0;
for (; isIt != isEndIt; ++ isIt) {
// store intersection, this might be costly
const auto& intersection = *isIt;
// deal with grid boundaries
if (intersection.boundary()) {
// compute the transmissibilty for the boundary intersection
const auto& geometry = intersection.geometry();
const auto& faceCenterInside = geometry.center();
auto faceAreaNormal = intersection.centerUnitOuterNormal();
faceAreaNormal *= geometry.volume();
Scalar transBoundaryIs;
computeHalfTrans_(transBoundaryIs,
faceAreaNormal,
intersection.indexInInside(),
distanceVector_(faceCenterInside,
intersection.indexInInside(),
elemIdx,
axisCentroids),
permeability_[elemIdx]);
// normally there would be two half-transmissibilities that would be
// averaged. on the grid boundary there only is the half
// transmissibility of the interior element.
transBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = transBoundaryIs;
// for boundary intersections we also need to compute the thermal
// half transmissibilities
if (enableEnergy) {
const auto& n = intersection.centerUnitOuterNormal();
const auto& inPos = elem.geometry().center();
const auto& outPos = intersection.geometry().center();
const auto& d = outPos - inPos;
// eWoms expects fluxes to be area specific, i.e. we must *not*
// the transmissibility with the face area here
Scalar thermalHalfTrans = std::abs(n*d)/(d*d);
thermalHalfTransBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] =
thermalHalfTrans;
}
++ boundaryIsIdx;
continue;
}
if (!intersection.neighbor())
// elements can be on process boundaries, i.e. they are not on the
// domain boundary yet they don't have neighbors.
continue;
const auto& outsideElem = intersection.outside();
unsigned outsideElemIdx = elemMapper.index(outsideElem);
// update the "thermal half transmissibility" for the intersection
if (enableEnergy) {
const auto& n = intersection.centerUnitOuterNormal();
Scalar A = intersection.geometry().volume();
const auto& inPos = elem.geometry().center();
const auto& outPos = intersection.geometry().center();
const auto& d = outPos - inPos;
(*thermalHalfTrans_)[directionalIsId_(elemIdx, outsideElemIdx)] =
A * (n*d)/(d*d);
}
// we only need to calculate a face's transmissibility
// once...
if (elemIdx > outsideElemIdx)
continue;
unsigned insideCartElemIdx = cartMapper.cartesianIndex(elemIdx);
unsigned outsideCartElemIdx = cartMapper.cartesianIndex(outsideElemIdx);
// local indices of the faces of the inside and
// outside elements which contain the intersection
int insideFaceIdx = intersection.indexInInside();
int outsideFaceIdx = intersection.indexInOutside();
if (insideFaceIdx == -1) {
// NNC. Set zero transmissibility, as it will be
// *added to* by applyNncToGridTrans_() later.
assert(outsideFaceIdx == -1);
trans_[isId_(elemIdx, outsideElemIdx)] = 0.0;
continue;
}
DimVector faceCenterInside;
DimVector faceCenterOutside;
DimVector faceAreaNormal;
typename std::is_same<Grid, Dune::CpGrid>::type isCpGrid;
computeFaceProperties(intersection,
elemIdx,
insideFaceIdx,
outsideElemIdx,
outsideFaceIdx,
faceCenterInside,
faceCenterOutside,
faceAreaNormal,
isCpGrid);
Scalar halfTrans1;
Scalar halfTrans2;
computeHalfTrans_(halfTrans1,
faceAreaNormal,
insideFaceIdx,
distanceVector_(faceCenterInside,
intersection.indexInInside(),
elemIdx,
axisCentroids),
permeability_[elemIdx]);
computeHalfTrans_(halfTrans2,
faceAreaNormal,
outsideFaceIdx,
distanceVector_(faceCenterOutside,
intersection.indexInOutside(),
outsideElemIdx,
axisCentroids),
permeability_[outsideElemIdx]);
applyNtg_(halfTrans1, insideFaceIdx, insideCartElemIdx, ntg);
applyNtg_(halfTrans2, outsideFaceIdx, outsideCartElemIdx, ntg);
// convert half transmissibilities to full face
// transmissibilities using the harmonic mean
Scalar trans;
if (std::abs(halfTrans1) < 1e-30 || std::abs(halfTrans2) < 1e-30)
// avoid division by zero
trans = 0.0;
else
trans = 1.0 / (1.0/halfTrans1 + 1.0/halfTrans2);
// apply the full face transmissibility multipliers
// for the inside ...
// The MULTZ needs special case if the option is ALL
// Then the smallest multiplier is applied.
// Default is to apply the top and bottom multiplier
bool useSmallestMultiplier = eclGrid.getMultzOption() == Opm::PinchMode::ModeEnum::ALL;
if (useSmallestMultiplier)
applyAllZMultipliers_(trans, insideFaceIdx, insideCartElemIdx, outsideCartElemIdx, transMult, cartDims);
else
applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
// ... and outside elements
applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult);
// apply the region multipliers (cf. the MULTREGT keyword)
Opm::FaceDir::DirEnum faceDir;
switch (insideFaceIdx) {
case 0:
case 1:
faceDir = Opm::FaceDir::XPlus;
break;
case 2:
case 3:
faceDir = Opm::FaceDir::YPlus;
break;
case 4:
case 5:
faceDir = Opm::FaceDir::ZPlus;
break;
default:
throw std::logic_error("Could not determine a face direction");
}
trans *= transMult.getRegionMultiplier(insideCartElemIdx,
outsideCartElemIdx,
faceDir);
trans_[isId_(elemIdx, outsideElemIdx)] = trans;
}
}
// potentially overwrite and/or modify transmissibilities based on input from deck
updateFromEclState_();
// Create mapping from global to local index
const size_t cartesianSize = cartMapper.cartesianSize();
// reserve memory
std::vector<int> globalToLocal(cartesianSize, -1);
// loop over all elements (global grid) and store Cartesian index
elemIt = vanguard_.grid().leafGridView().template begin<0>();
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = elemMapper.index(*elemIt);
int cartElemIdx = vanguard_.cartesianIndexMapper().cartesianIndex(elemIdx);
globalToLocal[cartElemIdx] = elemIdx;
}
applyEditNncToGridTrans_(elemMapper, globalToLocal);
applyNncToGridTrans_(globalToLocal);
//remove very small non-neighbouring transmissibilities
removeSmallNonCartesianTransmissibilities_();
}
/*!
* \brief Return the permeability for an element.
*/
const DimMatrix& permeability(unsigned elemIdx) const
{ return permeability_[elemIdx]; }
/*!
* \brief Return the transmissibility for the intersection between two elements.
*/
Scalar transmissibility(unsigned elemIdx1, unsigned elemIdx2) const
{ return trans_.at(isId_(elemIdx1, elemIdx2)); }
/*!
* \brief Return the transmissibility for a given boundary segment.
*/
Scalar transmissibilityBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
{ return transBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx)); }
/*!
* \brief Return the thermal "half transmissibility" for the intersection between two
* elements.
*
* The "half transmissibility" features all sub-expressions of the "thermal
* transmissibility" which can be precomputed, i.e. they are not dependent on the
* current solution:
*
* H_t = A * (n*d)/(d*d);
*
* where A is the area of the intersection between the inside and outside elements, n
* is the outer unit normal, and d is the distance between the center of the inside
* cell and the center of the intersection.
*/
Scalar thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const
{ return thermalHalfTrans_->at(directionalIsId_(insideElemIdx, outsideElemIdx)); }
Scalar thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const
{ return thermalHalfTransBoundary_.at(std::make_pair(insideElemIdx, boundaryFaceIdx)); }
private:
void removeSmallNonCartesianTransmissibilities_()
{
const auto& cartMapper = vanguard_.cartesianIndexMapper();
const auto& cartDims = cartMapper.cartesianDimensions();
for (auto&& trans: trans_) {
if (trans.second < transmissibilityThreshold_) {
const auto& id = trans.first;
const auto& elements = isIdReverse_(id);
int gc1 = std::min(cartMapper.cartesianIndex(elements.first), cartMapper.cartesianIndex(elements.second));
int gc2 = std::max(cartMapper.cartesianIndex(elements.first), cartMapper.cartesianIndex(elements.second));
// only adjust the NNCs
if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1])
continue;
//remove transmissibilities less than the threshold (by default 1e-6 in the deck's unit system)
trans.second = 0.0;
}
}
}
void applyAllZMultipliers_(Scalar& trans,
unsigned insideFaceIdx,
unsigned insideCartElemIdx,
unsigned outsideCartElemIdx,
const Opm::TransMult& transMult,
const std::array<int, dimWorld>& cartDims)
{
if (insideFaceIdx > 3) { // top or or bottom
Scalar mult = 1e20;
unsigned cartElemIdx = insideCartElemIdx;
// pick the smallest multiplier while looking down the pillar untill reaching the other end of the connection
// for the inbetween cells we apply it from both sides
while (cartElemIdx != outsideCartElemIdx) {
if (insideFaceIdx == 4 || cartElemIdx !=insideCartElemIdx )
mult = std::min(mult, transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZMinus));
if (insideFaceIdx == 5 || cartElemIdx !=insideCartElemIdx)
mult = std::min(mult, transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZPlus));
cartElemIdx += cartDims[0]*cartDims[1];
}
trans *= mult;
}
else
applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult);
}
void updateFromEclState_()
{
const auto& gridView = vanguard_.gridView();
const auto& cartMapper = vanguard_.cartesianIndexMapper();
const auto& cartDims = cartMapper.cartesianDimensions();
const auto& properties = vanguard_.eclState().get3DProperties();
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
#else
ElementMapper elemMapper(gridView);
#endif
const auto& inputTranx = properties.getDoubleGridProperty("TRANX");
const auto& inputTrany = properties.getDoubleGridProperty("TRANY");
const auto& inputTranz = properties.getDoubleGridProperty("TRANZ");
// compute the transmissibilities for all intersections
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
auto isIt = gridView.ibegin(elem);
const auto& isEndIt = gridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
// store intersection, this might be costly
const auto& intersection = *isIt;
if (!intersection.neighbor())
continue; // intersection is on the domain boundary
unsigned c1 = elemMapper.index(intersection.inside());
unsigned c2 = elemMapper.index(intersection.outside());
if (c1 > c2)
continue; // we only need to handle each connection once, thank you.
auto isId = isId_(c1, c2);
int gc1 = std::min(cartMapper.cartesianIndex(c1), cartMapper.cartesianIndex(c2));
int gc2 = std::max(cartMapper.cartesianIndex(c1), cartMapper.cartesianIndex(c2));
if (gc2 - gc1 == 1) {
if (inputTranx.deckAssigned())
// set simulator internal transmissibilities to values from inputTranx
trans_[isId] = inputTranx.iget(gc1);
else
// Scale transmissibilities with scale factor from inputTranx
trans_[isId] *= inputTranx.iget(gc1);
}
else if (gc2 - gc1 == cartDims[0]) {
if (inputTrany.deckAssigned())
// set simulator internal transmissibilities to values from inputTrany
trans_[isId] = inputTrany.iget(gc1);
else
// Scale transmissibilities with scale factor from inputTrany
trans_[isId] *= inputTrany.iget(gc1);
}
else if (gc2 - gc1 == cartDims[0]*cartDims[1]) {
if (inputTranz.deckAssigned())
// set simulator internal transmissibilities to values from inputTranz
trans_[isId] = inputTranz.iget(gc1);
else
// Scale transmissibilities with scale factor from inputTranz
trans_[isId] *= inputTranz.iget(gc1);
}
//else.. We don't support modification of NNC at the moment.
}
}
}
template <class Intersection>
void computeFaceProperties(const Intersection& intersection,
const int insideElemIdx,
const int insideFaceIdx,
const int outsideElemIdx,
const int outsideFaceIdx,
DimVector& faceCenterInside,
DimVector& faceCenterOutside,
DimVector& faceAreaNormal,
/*isCpGrid=*/std::false_type) const
{
// default implementation for DUNE grids
const auto& geometry = intersection.geometry();
faceCenterInside = geometry.center();
faceCenterOutside = faceCenterInside;
faceAreaNormal = intersection.centerUnitOuterNormal();
faceAreaNormal *= geometry.volume();
}
template <class Intersection>
void computeFaceProperties(const Intersection& intersection,
const int insideElemIdx,
const int insideFaceIdx,
const int outsideElemIdx,
const int outsideFaceIdx,
DimVector& faceCenterInside,
DimVector& faceCenterOutside,
DimVector& faceAreaNormal,
/*isCpGrid=*/std::true_type) const
{
int faceIdx = intersection.id();
faceCenterInside = vanguard_.grid().faceCenterEcl(insideElemIdx, insideFaceIdx);
faceCenterOutside = vanguard_.grid().faceCenterEcl(outsideElemIdx, outsideFaceIdx);
faceAreaNormal = vanguard_.grid().faceAreaNormalEcl(faceIdx);
}
/*
* \brief Applies additional transmissibilities specified via NNC keyword.
*
* Applies only those NNC that are actually represented by the grid. These may
* NNCs due to faults or NNCs that are actually neighbours. In both cases that
* specified transmissibilities (scaled by EDITNNC) will be added to the already
* existing models.
*
* \param cartesianToCompressed Vector containing the compressed index (or -1 for inactive
* cells) as the element at the cartesian index.
* \return Two vector of NNCs (scaled by EDITNNC). The first one are the NNCs that have been applied
* and the second the NNCs not resembled by faces of the grid. NNCs specified for
* inactive cells are omitted in these vectors.
*/
std::tuple<std::vector<Opm::NNCdata>, std::vector<Opm::NNCdata> >
applyNncToGridTrans_(const std::vector<int>& cartesianToCompressed)
{
// First scale NNCs with EDITNNC.
std::vector<Opm::NNCdata> unprocessedNnc;
std::vector<Opm::NNCdata> processedNnc;
const auto& nnc = vanguard_.eclState().getInputNNC();
if (!nnc.hasNNC())
return make_tuple(processedNnc, unprocessedNnc);
auto nncData = nnc.nncdata();
auto editnncData = vanguard_.eclState().getInputEDITNNC().data();
auto nncLess =
[](const Opm::NNCdata& d1, const Opm::NNCdata& d2)
{
return
(d1.cell1 < d2.cell1)
|| (d1.cell1 == d2.cell1 && d1.cell2 < d2.cell2);
};
std::sort(nncData.begin(), nncData.end(), nncLess);
auto candidate = nncData.begin();
for (const auto& edit: editnncData) {
auto printNncWarning =
[](int c1, int c2)
{
std::ostringstream sstr;
sstr << "Cannot edit NNC from " << c1 << " to " << c2
<< " as it does not exist";
Opm::OpmLog::warning(sstr.str());
};
if (candidate == nncData.end()) {
// no more NNCs left
printNncWarning(edit.cell1, edit.cell2);
continue;
}
if (candidate->cell1 != edit.cell1 || candidate->cell2 != edit.cell2) {
candidate = std::lower_bound(candidate, nncData.end(), Opm::NNCdata(edit.cell1, edit.cell2, 0), nncLess);
if (candidate == nncData.end()) {
// no more NNCs left
printNncWarning(edit.cell1, edit.cell2);
continue;
}
}
auto firstCandidate = candidate;
while (candidate != nncData.end()
&& candidate->cell1 == edit.cell1
&& candidate->cell2 == edit.cell2)
{
candidate->trans *= edit.trans;
++candidate;
}
// start with first match in next iteration to catch case where next
// EDITNNC is for same pair.
candidate = firstCandidate;
}
for (const auto& nncEntry : nnc.nncdata()) {
auto c1 = nncEntry.cell1;
auto c2 = nncEntry.cell2;
auto low = cartesianToCompressed[c1];
auto high = cartesianToCompressed[c2];
if (low > high)
std::swap(low, high);
if (low == -1 && high == -1)
// Silently discard as it is not between active cells
continue;
if (low == -1 || high == -1)
OPM_THROW(std::logic_error, "NNC between active and inactive cells (" <<
low << " -> " << high);
auto candidate = trans_.find(isId_(low, high));
if (candidate == trans_.end())
// This NNC is not resembled by the grid. Save it for later
// processing with local cell values
unprocessedNnc.push_back({c1, c2, nncEntry.trans});
else {
// NNC is represented by the grid and might be a neighboring connection
// In this case the transmissibilty is added to the value already
// set or computed.
candidate->second += nncEntry.trans;
processedNnc.push_back({c1, c2, nncEntry.trans});
}
}
return make_tuple(processedNnc, unprocessedNnc);
}
/// \brief Multiplies the grid transmissibilities according to EDITNNC.
void applyEditNncToGridTrans_(const ElementMapper& elementMapper,
const std::vector<int>& globalToLocal)
{
const auto& editNnc = vanguard_.eclState().getInputEDITNNC();
if (editNnc.empty())
return;
// editNnc is supposed to only reference non-neighboring connections and not
// neighboring connections. Use all entries for scaling if there is an NNC.
// variable nnc incremented in loop body.
auto nnc = editNnc.data().begin();
auto end = editNnc.data().end();
while (nnc != end) {
auto c1 = nnc->cell1;
auto c2 = nnc->cell2;
auto low = globalToLocal[c1];
auto high = globalToLocal[c2];
if (low > high)
std::swap(low, high);
auto candidate = trans_.find(isId_(low, high));
if (candidate == trans_.end()) {
std::ostringstream sstr;
sstr << "Cannot edit NNC from " << c1 << " to " << c2
<< " as it does not exist";
Opm::OpmLog::warning(sstr.str());
++nnc;
}
else {
// NNC exists
while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) {
candidate->second *= nnc->trans;
++nnc;
}
}
}
}
void extractPermeability_()
{
const auto& props = vanguard_.eclState().get3DProperties();
unsigned numElem = vanguard_.gridView().size(/*codim=*/0);
permeability_.resize(numElem);
// read the intrinsic permeabilities from the eclState. Note that all arrays
// provided by eclState are one-per-cell of "uncompressed" grid, whereas the
// simulation grid might remove a few elements. (e.g. because it is distributed
// over several processes.)
if (props.hasDeckDoubleGridProperty("PERMX")) {
const std::vector<double>& permxData =
props.getDoubleGridProperty("PERMX").getData();
std::vector<double> permyData(permxData);
if (props.hasDeckDoubleGridProperty("PERMY"))
permyData = props.getDoubleGridProperty("PERMY").getData();
std::vector<double> permzData(permxData);
if (props.hasDeckDoubleGridProperty("PERMZ"))
permzData = props.getDoubleGridProperty("PERMZ").getData();
for (size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
unsigned cartesianElemIdx = vanguard_.cartesianIndex(dofIdx);
permeability_[dofIdx] = 0.0;
permeability_[dofIdx][0][0] = permxData[cartesianElemIdx];
permeability_[dofIdx][1][1] = permyData[cartesianElemIdx];
permeability_[dofIdx][2][2] = permzData[cartesianElemIdx];
}
// for now we don't care about non-diagonal entries
}
else
throw std::logic_error("Can't read the intrinsic permeability from the ecl state. "
"(The PERM{X,Y,Z} keywords are missing)");
}
std::uint64_t isId_(std::uint32_t elemIdx1, std::uint32_t elemIdx2) const
{
std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2);
std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2);
return (elemBIdx<<elemIdxShift) + elemAIdx;
}
std::pair<std::uint32_t, std::uint32_t> isIdReverse_(const std::uint64_t& id) const
{
// Assigning an unsigned integer to a narrower type discards the most significant bits.
// See "The C programming language", section A.6.2.
// NOTE that the ordering of element A and B may have changed
std::uint32_t elemAIdx = id;
std::uint32_t elemBIdx = (id - elemAIdx) >> elemIdxShift;
return std::make_pair(elemAIdx, elemBIdx);
}
std::uint64_t directionalIsId_(std::uint32_t elemIdx1, std::uint32_t elemIdx2) const
{
return (std::uint64_t(elemIdx1)<<elemIdxShift) + elemIdx2;
}
void computeHalfTrans_(Scalar& halfTrans,
const DimVector& areaNormal,
int faceIdx, // in the reference element that contains the intersection
const DimVector& distance,
const DimMatrix& perm) const
{
assert(faceIdx >= 0);
unsigned dimIdx = faceIdx/2;
assert(dimIdx < dimWorld);
halfTrans = perm[dimIdx][dimIdx];
Scalar val = 0;
for (unsigned i = 0; i < areaNormal.size(); ++i)
val += areaNormal[i]*distance[i];
halfTrans *= std::abs(val);
halfTrans /= distance.two_norm2();
}
DimVector distanceVector_(const DimVector& center,
int faceIdx, // in the reference element that contains the intersection
unsigned elemIdx,
const std::array<std::vector<DimVector>, dimWorld>& axisCentroids) const
{
assert(faceIdx >= 0);
unsigned dimIdx = faceIdx/2;
assert(dimIdx < dimWorld);
DimVector x = center;
x -= axisCentroids[dimIdx][elemIdx];
return x;
}
void applyMultipliers_(Scalar& trans,
unsigned faceIdx,
unsigned cartElemIdx,
const Opm::TransMult& transMult) const
{
// apply multiplyer for the transmissibility of the face. (the
// face index is the index of the reference-element face which
// contains the intersection of interest.)
switch (faceIdx) {
case 0: // left
trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::XMinus);
break;
case 1: // right
trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::XPlus);
break;
case 2: // front
trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::YMinus);
break;
case 3: // back
trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::YPlus);
break;
case 4: // bottom
trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZMinus);
break;
case 5: // top
trans *= transMult.getMultiplier(cartElemIdx, Opm::FaceDir::ZPlus);
break;
}
}
void applyNtg_(Scalar& trans,
unsigned faceIdx,
unsigned cartElemIdx,
const std::vector<double>& ntg) const
{
// apply multiplyer for the transmissibility of the face. (the
// face index is the index of the reference-element face which
// contains the intersection of interest.)
switch (faceIdx) {
case 0: // left
trans *= ntg[cartElemIdx];
break;
case 1: // right
trans *= ntg[cartElemIdx];
break;
case 2: // front
trans *= ntg[cartElemIdx];
break;
case 3: // back
trans *= ntg[cartElemIdx];
break;
// NTG does not apply to top and bottom faces
}
}
void minPvFillNtg_(std::vector<double>& averageNtg) const
{
// compute volume weighted arithmetic average of NTG for
// cells merged as an result of minpv.
const auto& eclState = vanguard_.eclState();
const auto& eclGrid = eclState.getInputGrid();
const std::vector<double>& ntg =
eclState.get3DProperties().getDoubleGridProperty("NTG").getData();
averageNtg = ntg;
bool opmfil = eclGrid.getMinpvMode() == Opm::MinpvMode::OpmFIL;
// just return the unmodified ntg if opmfil is not used
if (!opmfil)
return;
const auto& porv = eclState.get3DProperties().getDoubleGridProperty("PORV").getData();
const auto& actnum = eclState.get3DProperties().getIntGridProperty("ACTNUM").getData();
const auto& cartMapper = vanguard_.cartesianIndexMapper();
const auto& cartDims = cartMapper.cartesianDimensions();
assert(dimWorld > 1);
const size_t nxny = cartDims[0] * cartDims[1];
for (size_t cartesianCellIdx = 0; cartesianCellIdx < ntg.size(); ++cartesianCellIdx) {
// use the original ntg values for the inactive cells
if (!actnum[cartesianCellIdx])
continue;
// Average properties as long as there exist cells above
// that has pore volume less than the MINPV threshold
const double cellVolume = eclGrid.getCellVolume(cartesianCellIdx);
double ntgCellVolume = ntg[cartesianCellIdx] * cellVolume;
double totalCellVolume = cellVolume;
int cartesianCellIdxAbove = cartesianCellIdx - nxny;
while (cartesianCellIdxAbove >= 0 &&
actnum[cartesianCellIdxAbove] > 0 &&
porv[cartesianCellIdxAbove] < eclGrid.getMinpvVector()[cartesianCellIdxAbove])
{
// Volume weighted arithmetic average of NTG
const double cellAboveVolume = eclGrid.getCellVolume(cartesianCellIdxAbove);
totalCellVolume += cellAboveVolume;
ntgCellVolume += ntg[cartesianCellIdxAbove]*cellAboveVolume;
cartesianCellIdxAbove -= nxny;
}
averageNtg[cartesianCellIdx] = ntgCellVolume / totalCellVolume;
}
}
const Vanguard& vanguard_;
Scalar transmissibilityThreshold_;
std::vector<DimMatrix> permeability_;
std::unordered_map<std::uint64_t, Scalar> trans_;
std::map<std::pair<unsigned, unsigned>, Scalar> transBoundary_;
std::map<std::pair<unsigned, unsigned>, Scalar> thermalHalfTransBoundary_;
Opm::ConditionalStorage<enableEnergy,
std::unordered_map<std::uint64_t, Scalar> > thermalHalfTrans_;
};
} // namespace Ewoms
#endif

821
ebos/eclwellmanager.hh Normal file
View File

@ -0,0 +1,821 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/**
* \file
*
* \copydoc Ewoms::EclWellManager
*/
#ifndef EWOMS_ECL_WELL_MANAGER_HH
#define EWOMS_ECL_WELL_MANAGER_HH
#include "eclpeacemanwell.hh"
#include <ewoms/disc/common/fvbaseproperties.hh>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellConnections.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/TimeMap.hpp>
#include <opm/output/data/Wells.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <ewoms/common/propertysystem.hh>
#include <ewoms/parallel/threadedentityiterator.hh>
#include <dune/grid/common/gridenums.hh>
#include <map>
#include <string>
#include <vector>
BEGIN_PROPERTIES
NEW_PROP_TAG(Grid);
END_PROPERTIES
namespace Ewoms {
/*!
* \ingroup EclBlackOilSimulator
*
* \brief A class which handles well controls as specified by an
* Eclipse deck
*/
template <class TypeTag>
class EclWellManager
{
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
enum { numPhases = FluidSystem::numPhases };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
typedef typename GridView::template Codim<0>::Entity Element;
typedef Ewoms::EclPeacemanWell<TypeTag> Well;
typedef std::map<int, std::pair<const Opm::Connection*, std::shared_ptr<Well> > > WellConnectionsMap;
typedef Dune::FieldVector<Evaluation, numEq> EvalEqVector;
public:
EclWellManager(Simulator& simulator)
: simulator_(simulator)
{ }
/*!
* \brief This sets up the basic properties of all wells.
*
* I.e., well positions, names etc...
*/
void init(const Opm::EclipseState& eclState OPM_UNUSED,
const Opm::Schedule& deckSchedule)
{
// create the wells which intersect with the current process' grid
for (size_t deckWellIdx = 0; deckWellIdx < deckSchedule.numWells(); ++deckWellIdx)
{
const Opm::Well* deckWell = deckSchedule.getWells()[deckWellIdx];
const std::string& wellName = deckWell->name();
Scalar wellTemperature = 273.15 + 15.56; // [K]
if (deckWell->isInjector(/*timeStep=*/0))
wellTemperature = deckWell->getInjectionProperties(/*timeStep=*/0).temperature;
// set the name of the well but not much else. (i.e., if it is not completed,
// the well primarily serves as a placeholder.) The big rest of the well is
// specified by the updateWellCompletions_() method
auto well = std::make_shared<Well>(simulator_);
well->setName(wellName);
well->setWellStatus(Well::Shut);
well->setTemperature(wellTemperature);
wells_.push_back(well);
wellNameToIndex_[well->name()] = wells_.size() - 1;
}
}
/*!
* \brief This should be called the problem before each simulation
* episode to adapt the well controls.
*/
void beginEpisode(const Opm::EclipseState& eclState, const Opm::Schedule& deckSchedule, bool wasRestarted=false)
{
unsigned episodeIdx = simulator_.episodeIndex();
WellConnectionsMap wellCompMap;
computeWellConnectionsMap_(episodeIdx, wellCompMap);
if (wasRestarted || wellTopologyChanged_(eclState, deckSchedule, episodeIdx))
updateWellTopology_(episodeIdx, wellCompMap, gridDofIsPenetrated_);
// set those parameters of the wells which do not change the topology of the
// linearized system of equations
updateWellParameters_(episodeIdx, wellCompMap);
const std::vector<const Opm::Well*>& deckWells = deckSchedule.getWells(episodeIdx);
// set the injection data for the respective wells.
for (size_t deckWellIdx = 0; deckWellIdx < deckWells.size(); ++deckWellIdx) {
const Opm::Well* deckWell = deckWells[deckWellIdx];
if (!hasWell(deckWell->name()))
continue;
auto well = this->well(deckWell->name());
if (deckWell->isInjector(episodeIdx))
well->setTemperature(deckWell->getInjectionProperties(episodeIdx).temperature);
Opm::WellCommon::StatusEnum deckWellStatus = deckWell->getStatus(episodeIdx);
switch (deckWellStatus) {
case Opm::WellCommon::AUTO:
// TODO: for now, auto means open...
case Opm::WellCommon::OPEN:
well->setWellStatus(Well::Open);
break;
case Opm::WellCommon::STOP:
well->setWellStatus(Well::Closed);
break;
case Opm::WellCommon::SHUT:
well->setWellStatus(Well::Shut);
break;
}
// make sure that the well is either an injector or a
// producer for the current episode. (it is not allowed to
// be neither or to be both...)
assert((deckWell->isInjector(episodeIdx)?1:0) +
(deckWell->isProducer(episodeIdx)?1:0) == 1);
if (deckWell->isInjector(episodeIdx)) {
well->setWellType(Well::Injector);
const Opm::WellInjectionProperties& injectProperties =
deckWell->getInjectionProperties(episodeIdx);
switch (injectProperties.injectorType) {
case Opm::WellInjector::WATER:
well->setInjectedPhaseIndex(FluidSystem::waterPhaseIdx);
break;
case Opm::WellInjector::GAS:
well->setInjectedPhaseIndex(FluidSystem::gasPhaseIdx);
break;
case Opm::WellInjector::OIL:
well->setInjectedPhaseIndex(FluidSystem::oilPhaseIdx);
break;
case Opm::WellInjector::MULTI:
throw std::runtime_error("Not implemented: Multi-phase injector wells");
}
switch (injectProperties.controlMode) {
case Opm::WellInjector::RATE:
well->setControlMode(Well::ControlMode::VolumetricSurfaceRate);
break;
case Opm::WellInjector::RESV:
well->setControlMode(Well::ControlMode::VolumetricReservoirRate);
break;
case Opm::WellInjector::BHP:
well->setControlMode(Well::ControlMode::BottomHolePressure);
break;
case Opm::WellInjector::THP:
well->setControlMode(Well::ControlMode::TubingHeadPressure);
break;
case Opm::WellInjector::GRUP:
throw std::runtime_error("Not implemented: Well groups");
case Opm::WellInjector::CMODE_UNDEFINED:
std::cout << "Warning: Control mode of injection well " << well->name()
<< " is undefined. Assuming well to be shut.\n";
well->setWellStatus(Well::WellStatus::Shut);
continue;
}
switch (injectProperties.injectorType) {
case Opm::WellInjector::WATER:
well->setVolumetricPhaseWeights(/*oil=*/0.0, /*gas=*/0.0, /*water=*/1.0);
break;
case Opm::WellInjector::OIL:
well->setVolumetricPhaseWeights(/*oil=*/1.0, /*gas=*/0.0, /*water=*/0.0);
break;
case Opm::WellInjector::GAS:
well->setVolumetricPhaseWeights(/*oil=*/0.0, /*gas=*/1.0, /*water=*/0.0);
break;
case Opm::WellInjector::MULTI:
throw std::runtime_error("Not implemented: Multi-phase injection wells");
}
well->setMaximumSurfaceRate(injectProperties.surfaceInjectionRate);
well->setMaximumReservoirRate(injectProperties.reservoirInjectionRate);
well->setTargetBottomHolePressure(injectProperties.BHPLimit);
// TODO
well->setTargetTubingHeadPressure(1e30);
//well->setTargetTubingHeadPressure(injectProperties.THPLimit);
}
if (deckWell->isProducer(episodeIdx)) {
well->setWellType(Well::Producer);
const Opm::WellProductionProperties& producerProperties =
deckWell->getProductionProperties(episodeIdx);
switch (producerProperties.controlMode) {
case Opm::WellProducer::ORAT:
well->setControlMode(Well::ControlMode::VolumetricSurfaceRate);
well->setVolumetricPhaseWeights(/*oil=*/1.0, /*gas=*/0.0, /*water=*/0.0);
well->setMaximumSurfaceRate(producerProperties.OilRate);
break;
case Opm::WellProducer::GRAT:
well->setControlMode(Well::ControlMode::VolumetricSurfaceRate);
well->setVolumetricPhaseWeights(/*oil=*/0.0, /*gas=*/1.0, /*water=*/0.0);
well->setMaximumSurfaceRate(producerProperties.GasRate);
break;
case Opm::WellProducer::WRAT:
well->setControlMode(Well::ControlMode::VolumetricSurfaceRate);
well->setVolumetricPhaseWeights(/*oil=*/0.0, /*gas=*/0.0, /*water=*/1.0);
well->setMaximumSurfaceRate(producerProperties.WaterRate);
break;
case Opm::WellProducer::LRAT:
well->setControlMode(Well::ControlMode::VolumetricSurfaceRate);
well->setVolumetricPhaseWeights(/*oil=*/1.0, /*gas=*/0.0, /*water=*/1.0);
well->setMaximumSurfaceRate(producerProperties.LiquidRate);
break;
case Opm::WellProducer::CRAT:
throw std::runtime_error("Not implemented: Linearly combined rates");
case Opm::WellProducer::RESV:
well->setControlMode(Well::ControlMode::VolumetricReservoirRate);
well->setVolumetricPhaseWeights(/*oil=*/1.0, /*gas=*/1.0, /*water=*/1.0);
well->setMaximumSurfaceRate(producerProperties.ResVRate);
break;
case Opm::WellProducer::BHP:
well->setControlMode(Well::ControlMode::BottomHolePressure);
break;
case Opm::WellProducer::THP:
well->setControlMode(Well::ControlMode::TubingHeadPressure);
break;
case Opm::WellProducer::GRUP:
throw std::runtime_error("Not implemented: Well groups");
case Opm::WellProducer::NONE:
// fall-through
case Opm::WellProducer::CMODE_UNDEFINED:
std::cout << "Warning: Control mode of production well " << well->name()
<< " is undefined. Assuming well to be shut.";
well->setWellStatus(Well::WellStatus::Shut);
continue;
}
well->setTargetBottomHolePressure(producerProperties.BHPLimit);
// TODO
well->setTargetTubingHeadPressure(-1e30);
//well->setTargetTubingHeadPressure(producerProperties.THPLimit);
}
}
}
/*!
* \brief Return the number of wells considered by the EclWellManager.
*/
unsigned numWells() const
{ return wells_.size(); }
/*!
* \brief Return if a given well name is known to the wells manager
*/
bool hasWell(const std::string& wellName) const
{
return wellNameToIndex_.find( wellName ) != wellNameToIndex_.end();
}
/*!
* \brief Returns true iff a given degree of freedom is currently penetrated by any well.
*/
bool gridDofIsPenetrated(unsigned globalDofIdx) const
{ return gridDofIsPenetrated_[globalDofIdx]; }
/*!
* \brief Given a well name, return the corresponding index.
*
* A std::runtime_error will be thrown if the well name is unknown.
*/
unsigned wellIndex(const std::string& wellName) const
{
assert( hasWell( wellName ) );
const auto& it = wellNameToIndex_.find(wellName);
if (it == wellNameToIndex_.end())
throw std::runtime_error("No well called '"+wellName+"'found");
return it->second;
}
/*!
* \brief Given a well name, return the corresponding well.
*
* A std::runtime_error will be thrown if the well name is unknown.
*/
std::shared_ptr<const Well> well(const std::string& wellName) const
{ return wells_[wellIndex(wellName)]; }
/*!
* \brief Given a well name, return the corresponding well.
*
* A std::runtime_error will be thrown if the well name is unknown.
*/
std::shared_ptr<Well> well(const std::string& wellName)
{ return wells_[wellIndex(wellName)]; }
/*!
* \brief Given a well index, return the corresponding well.
*/
std::shared_ptr<const Well> well(size_t wellIdx) const
{ return wells_[wellIdx]; }
/*!
* \brief Given a well index, return the corresponding well.
*/
std::shared_ptr<Well> well(size_t wellIdx)
{ return wells_[wellIdx]; }
/*!
* \brief Informs the well manager that a time step has just begun.
*/
void beginTimeStep()
{
// iterate over all wells and notify them individually
for (size_t wellIdx = 0; wellIdx < wells_.size(); ++wellIdx)
wells_[wellIdx]->beginTimeStep();
}
/*!
* \brief Informs the well that an iteration has just begun.
*
* In this method, the well calculates the bottom hole and tubing head pressures, the
* actual unconstraint production and injection rates, etc.
*/
void beginIteration()
{
// call the preprocessing routines
const size_t wellSize = wells_.size();
for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx)
wells_[wellIdx]->beginIterationPreProcess();
// call the accumulation routines
ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(simulator_.vanguard().gridView());
#ifdef _OPENMP
#pragma omp parallel
#endif
{
ElementContext elemCtx(simulator_);
auto elemIt = threadedElemIt.beginParallel();
for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
const Element& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity)
continue;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx)
wells_[wellIdx]->beginIterationAccumulate(elemCtx, /*timeIdx=*/0);
}
}
// call the postprocessing routines
for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx)
wells_[wellIdx]->beginIterationPostProcess();
}
/*!
* \brief Informs the well manager that an iteration has just been finished.
*/
void endIteration()
{
// iterate over all wells and notify them individually
const size_t wellSize = wells_.size();
for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx)
wells_[wellIdx]->endIteration();
}
/*!
* \brief Informs the well manager that a time step has just been finished.
*/
void endTimeStep()
{
Scalar dt = simulator_.timeStepSize();
// iterate over all wells and notify them individually. also, update the
// production/injection totals for the active wells.
const size_t wellSize = wells_.size();
for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx) {
auto well = wells_[wellIdx];
well->endTimeStep();
// update the surface volumes of the produced/injected fluids
std::array<Scalar, numPhases>* injectedVolume;
if (wellTotalInjectedVolume_.count(well->name()) == 0) {
injectedVolume = &wellTotalInjectedVolume_[well->name()];
std::fill(injectedVolume->begin(), injectedVolume->end(), 0.0);
}
else
injectedVolume = &wellTotalInjectedVolume_[well->name()];
std::array<Scalar, numPhases>* producedVolume;
if (wellTotalProducedVolume_.count(well->name()) == 0) {
producedVolume = &wellTotalProducedVolume_[well->name()];
std::fill(producedVolume->begin(), producedVolume->end(), 0.0);
}
else
producedVolume = &wellTotalProducedVolume_[well->name()];
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
// this assumes that the implicit Euler method is used for time
// integration. TODO: Once the time discretization becomes pluggable,
// this integration needs to be done by the time discretization code!
Scalar vol = dt * well->surfaceRate(phaseIdx);
if (vol < 0)
(*producedVolume)[phaseIdx] += -vol;
else
(*injectedVolume)[phaseIdx] += vol;
}
}
}
/*!
* \brief Informs the well manager that a simulation episode has just been finished.
*/
void endEpisode()
{ }
/*!
* \brief Returns the surface volume of a fluid phase produced by a well.
*/
Scalar totalProducedVolume(const std::string& wellName, unsigned phaseIdx) const
{
if (wellTotalProducedVolume_.count(wellName) == 0)
return 0.0; // well not yet seen
return wellTotalProducedVolume_.at(wellName)[phaseIdx];
}
/*!
* \brief Returns the surface volume of a fluid phase injected by a well.
*/
Scalar totalInjectedVolume(const std::string& wellName, unsigned phaseIdx) const
{
if (wellTotalInjectedVolume_.count(wellName) == 0)
return 0.0; // well not yet seen
return wellTotalInjectedVolume_.at(wellName)[phaseIdx];
}
/*!
* \brief Computes the source term due to wells for a degree of
* freedom.
*/
template <class Context>
void computeTotalRatesForDof(EvalEqVector& q,
const Context& context,
unsigned dofIdx,
unsigned timeIdx) const
{
q = 0.0;
if (!gridDofIsPenetrated(context.globalSpaceIndex(dofIdx, timeIdx)))
return;
RateVector wellRate;
// iterate over all wells and add up their individual rates
const size_t wellSize = wells_.size();
for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx) {
wellRate = 0.0;
wells_[wellIdx]->computeTotalRatesForDof(wellRate, context, dofIdx, timeIdx);
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
q[eqIdx] += wellRate[eqIdx];
}
}
Opm::data::Wells wellData() const
{
Opm::data::Wells wellDat;
using rt = Opm::data::Rates::opt;
for (unsigned wellIdx = 0; wellIdx < numWells(); ++wellIdx) {
const auto& ebosWell = well(wellIdx);
auto& wellOut = wellDat[ebosWell->name()];
wellOut.bhp = ebosWell->bottomHolePressure();
wellOut.thp = ebosWell->tubingHeadPressure();
wellOut.temperature = 0;
wellOut.rates.set( rt::wat, ebosWell->surfaceRate(waterPhaseIdx) );
wellOut.rates.set( rt::oil, ebosWell->surfaceRate(oilPhaseIdx) );
wellOut.rates.set( rt::gas, ebosWell->surfaceRate(gasPhaseIdx) );
const int numConnections = ebosWell->numConnections();
wellOut.connections.resize(numConnections);
for( int i = 0; i < numConnections; ++i ) {
auto& connection = wellOut.connections[ i ];
connection.index = 0;
connection.pressure = 0.0;
connection.reservoir_rate = 0.0;
connection.rates.set( rt::wat, 0.0 );
connection.rates.set( rt::oil, 0.0 );
connection.rates.set( rt::gas, 0.0 );
}
}
return wellDat;
}
/*!
* \brief This method writes the complete state of all wells
* to the hard disk.
*/
template <class Restarter>
void serialize(Restarter& res OPM_UNUSED)
{
/* do nothing: Everything which we need here is provided by the deck->.. */
}
/*!
* \brief This method restores the complete state of the all wells
* from disk.
*
* It is the inverse of the serialize() method.
*/
template <class Restarter>
void deserialize(Restarter& res OPM_UNUSED)
{
// initialize the wells for the current episode
beginEpisode(simulator_.vanguard().eclState(), simulator_.vanguard().schedule(), /*wasRestarted=*/true);
}
/*!
* \brief Returns true if something in a well changed compared to the previous report
* step.
*
* "Something" can either be the well topology (i.e., which grid blocks are contained
* in which well) or it can be a well parameter like the bottom hole pressure...
*/
bool wellsChanged(const Opm::EclipseState& eclState,
const Opm::Schedule& schedule,
unsigned reportStepIdx) const
{
if (wellTopologyChanged_(eclState, reportStepIdx))
return true;
if (schedule.getTimeMap().numTimesteps() <= (unsigned) reportStepIdx)
// for the "until the universe dies" episode, the wells don't change
return false;
const Opm::Events& events = schedule.getEvents();
return events.hasEvent(Opm::ScheduleEvents::PRODUCTION_UPDATE |
Opm::ScheduleEvents::INJECTION_UPDATE |
Opm::ScheduleEvents::WELL_STATUS_CHANGE,
reportStepIdx);
}
protected:
bool wellTopologyChanged_(const Opm::EclipseState& eclState OPM_UNUSED,
const Opm::Schedule& schedule,
unsigned reportStepIdx) const
{
if (reportStepIdx == 0) {
// the well topology has always been changed relative to before the
// simulation is started...
return true;
}
if (schedule.getTimeMap().numTimesteps() <= (unsigned) reportStepIdx)
// for the "until the universe dies" episode, the wells don't change
return false;
const Opm::Events& events = schedule.getEvents();
return events.hasEvent(Opm::ScheduleEvents::NEW_WELL |
Opm::ScheduleEvents::COMPLETION_CHANGE,
reportStepIdx);
}
void updateWellTopology_(unsigned reportStepIdx OPM_UNUSED,
const WellConnectionsMap& wellConnections,
std::vector<bool>& gridDofIsPenetrated) const
{
auto& model = simulator_.model();
const auto& vanguard = simulator_.vanguard();
// first, remove all wells from the reservoir
model.clearAuxiliaryModules();
auto wellIt = wells_.begin();
const auto& wellEndIt = wells_.end();
for (; wellIt != wellEndIt; ++wellIt) {
(*wellIt)->clear();
(*wellIt)->beginSpec();
}
//////
// tell the active wells which DOFs they contain
const auto gridView = simulator_.vanguard().gridView();
gridDofIsPenetrated.resize(model.numGridDof());
std::fill(gridDofIsPenetrated.begin(), gridDofIsPenetrated.end(), false);
ElementContext elemCtx(simulator_);
auto elemIt = gridView.template begin</*codim=*/0>();
const auto elemEndIt = gridView.template end</*codim=*/0>();
std::set<std::shared_ptr<Well> > wells;
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity)
continue; // non-local entities need to be skipped
elemCtx.updateStencil(elem);
for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx) {
unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
unsigned cartesianDofIdx = vanguard.cartesianIndex(globalDofIdx);
if (wellConnections.count(cartesianDofIdx) == 0)
// the current DOF is not contained in any well, so we must skip
// it...
continue;
gridDofIsPenetrated[globalDofIdx] = true;
auto eclWell = wellConnections.at(cartesianDofIdx).second;
eclWell->addDof(elemCtx, dofIdx);
wells.insert(eclWell);
}
//////
}
// register all wells at the model as auxiliary equations
wellIt = wells_.begin();
for (; wellIt != wellEndIt; ++wellIt) {
(*wellIt)->endSpec();
model.addAuxiliaryModule(wellIt->get());
}
}
void computeWellConnectionsMap_(unsigned reportStepIdx OPM_UNUSED, WellConnectionsMap& cartesianIdxToConnectionMap)
{
const auto& deckSchedule = simulator_.vanguard().schedule();
#ifndef NDEBUG
const auto& eclState = simulator_.vanguard().eclState();
const auto& eclGrid = eclState.getInputGrid();
assert( int(eclGrid.getNX()) == simulator_.vanguard().cartesianDimensions()[ 0 ] );
assert( int(eclGrid.getNY()) == simulator_.vanguard().cartesianDimensions()[ 1 ] );
assert( int(eclGrid.getNZ()) == simulator_.vanguard().cartesianDimensions()[ 2 ] );
#endif
// compute the mapping from logically Cartesian indices to the well the
// respective connection.
const std::vector<const Opm::Well*>& deckWells = deckSchedule.getWells(reportStepIdx);
for (size_t deckWellIdx = 0; deckWellIdx < deckWells.size(); ++deckWellIdx) {
const Opm::Well* deckWell = deckWells[deckWellIdx];
const std::string& wellName = deckWell->name();
if (!hasWell(wellName))
{
#ifndef NDEBUG
if( simulator_.vanguard().grid().comm().size() == 1 )
{
std::cout << "Well '" << wellName << "' suddenly appears in the connection "
<< "for the report step, but has not been previously specified. "
<< "Ignoring.\n";
}
#endif
continue;
}
std::array<int, 3> cartesianCoordinate;
// set the well parameters defined by the current set of connections
const auto& connectionSet = deckWell->getConnections(reportStepIdx);
for (size_t connIdx = 0; connIdx < connectionSet.size(); connIdx ++) {
const auto& connection = connectionSet.get(connIdx);
cartesianCoordinate[ 0 ] = connection.getI();
cartesianCoordinate[ 1 ] = connection.getJ();
cartesianCoordinate[ 2 ] = connection.getK();
unsigned cartIdx = simulator_.vanguard().cartesianIndex( cartesianCoordinate );
// in this code we only support each cell to be part of at most a single
// well. TODO (?) change this?
assert(cartesianIdxToConnectionMap.count(cartIdx) == 0);
auto eclWell = wells_[wellIndex(wellName)];
cartesianIdxToConnectionMap[cartIdx] = std::make_pair(&connection, eclWell);
}
}
}
void updateWellParameters_(unsigned reportStepIdx, const WellConnectionsMap& wellConnections)
{
const auto& deckSchedule = simulator_.vanguard().schedule();
const std::vector<const Opm::Well*>& deckWells = deckSchedule.getWells(reportStepIdx);
// set the reference depth for all wells
for (size_t deckWellIdx = 0; deckWellIdx < deckWells.size(); ++deckWellIdx) {
const Opm::Well* deckWell = deckWells[deckWellIdx];
const std::string& wellName = deckWell->name();
if( hasWell( wellName ) )
{
wells_[wellIndex(wellName)]->clear();
wells_[wellIndex(wellName)]->setReferenceDepth(deckWell->getRefDepth());
}
}
// associate the well connections with grid cells and register them in the
// Peaceman well object
const auto& vanguard = simulator_.vanguard();
const GridView gridView = vanguard.gridView();
ElementContext elemCtx(simulator_);
auto elemIt = gridView.template begin</*codim=*/0>();
const auto elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity)
continue; // non-local entities need to be skipped
elemCtx.updateStencil(elem);
for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx)
{
assert( elemCtx.numPrimaryDof(/*timeIdx=*/0) == 1 );
unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
unsigned cartesianDofIdx = vanguard.cartesianIndex(globalDofIdx);
if (wellConnections.count(cartesianDofIdx) == 0)
// the current DOF is not contained in any well, so we must skip
// it...
continue;
const auto& connInfo = wellConnections.at(cartesianDofIdx);
const Opm::Connection* connection = connInfo.first;
std::shared_ptr<Well> eclWell = connInfo.second;
eclWell->addDof(elemCtx, dofIdx);
eclWell->setConnectionTransmissibilityFactor(elemCtx, dofIdx, connection->CF());
eclWell->setRadius(elemCtx, dofIdx, connection->rw());
//eclWell->setEffectivePermeability(elemCtx, dofIdx, connection->Kh());
}
}
}
Simulator& simulator_;
std::vector<std::shared_ptr<Well> > wells_;
std::vector<bool> gridDofIsPenetrated_;
std::map<std::string, int> wellNameToIndex_;
std::map<std::string, std::array<Scalar, numPhases> > wellTotalInjectedVolume_;
std::map<std::string, std::array<Scalar, numPhases> > wellTotalProducedVolume_;
};
} // namespace Ewoms
#endif

497
ebos/eclwriter.hh Normal file
View File

@ -0,0 +1,497 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Ewoms::EclWriter
*/
#ifndef EWOMS_ECL_WRITER_HH
#define EWOMS_ECL_WRITER_HH
#include "collecttoiorank.hh"
#include "ecloutputblackoilmodule.hh"
#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
#include <ewoms/io/baseoutputwriter.hh>
#include <ewoms/parallel/tasklets.hh>
#include <opm/output/eclipse/EclipseIO.hpp>
#include <opm/output/eclipse/RestartValue.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <opm/grid/GridHelpers.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <list>
#include <utility>
#include <string>
BEGIN_PROPERTIES
NEW_PROP_TAG(EnableEclOutput);
NEW_PROP_TAG(EnableAsyncEclOutput);
NEW_PROP_TAG(EclOutputDoublePrecision);
END_PROPERTIES
namespace Ewoms {
template <class TypeTag>
class EclWriter;
template <class TypeTag>
class EclOutputBlackOilModule;
/*!
* \ingroup EclBlackOilSimulator
*
* \brief Collects necessary output values and pass it to opm-output.
*
* Caveats:
* - For this class to do do anything meaningful, you will have to
* have the OPM module opm-output.
* - The only DUNE grid which is currently supported is Dune::CpGrid
* from the OPM module "opm-grid". Using another grid won't
* fail at compile time but you will provoke a fatal exception as
* soon as you try to write an ECL output file.
* - This class requires to use the black oil model with the element
* centered finite volume discretization.
*/
template <class TypeTag>
class EclWriter
{
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef CollectDataToIORank<Vanguard> CollectDataToIORankType;
typedef std::vector<Scalar> ScalarBuffer;
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
public:
static void registerParameters()
{
EclOutputBlackOilModule<TypeTag>::registerParameters();
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableAsyncEclOutput,
"Write the ECL-formated results in a non-blocking way (i.e., using a separate thread).");
}
EclWriter(const Simulator& simulator)
: simulator_(simulator)
, collectToIORank_(simulator_.vanguard())
, eclOutputModule_(simulator, collectToIORank_)
{
globalGrid_ = simulator_.vanguard().grid();
globalGrid_.switchToGlobalView();
eclIO_.reset(new Opm::EclipseIO(simulator_.vanguard().eclState(),
Opm::UgGridHelpers::createEclipseGrid(globalGrid_, simulator_.vanguard().eclState().getInputGrid()),
simulator_.vanguard().schedule(),
simulator_.vanguard().summaryConfig()));
// create output thread if enabled and rank is I/O rank
// async output is enabled by default if pthread are enabled
bool enableAsyncOutput = EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput);
int numWorkerThreads = 0;
if (enableAsyncOutput && collectToIORank_.isIORank())
numWorkerThreads = 1;
taskletRunner_.reset(new TaskletRunner(numWorkerThreads));
}
~EclWriter()
{ }
const Opm::EclipseIO& eclIO() const
{ return *eclIO_; }
void writeInit()
{
if (collectToIORank_.isIORank()) {
std::map<std::string, std::vector<int> > integerVectors;
if (collectToIORank_.isParallel())
integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks());
eclIO_->writeInitial(computeTrans_(), integerVectors, exportNncStructure_());
}
}
/*!
* \brief collect and pass data and pass it to eclIO writer
*/
void writeOutput(bool isSubStep)
{
Scalar curTime = simulator_.time() + simulator_.timeStepSize();
Scalar totalCpuTime =
simulator_.executionTimer().realTimeElapsed() +
simulator_.setupTimer().realTimeElapsed() +
simulator_.vanguard().externalSetupTime();
Scalar nextStepSize = simulator_.problem().nextTimeStepSize();
// output using eclWriter if enabled
Opm::data::Wells localWellData = simulator_.problem().wellModel().wellData();
int episodeIdx = simulator_.episodeIndex() + 1;
const auto& gridView = simulator_.vanguard().gridView();
int numElements = gridView.size(/*codim=*/0);
bool log = collectToIORank_.isIORank();
eclOutputModule_.allocBuffers(numElements, episodeIdx, isSubStep, log);
ElementContext elemCtx(simulator_);
ElementIterator elemIt = gridView.template begin</*codim=*/0>();
const ElementIterator& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
eclOutputModule_.processElement(elemCtx);
}
eclOutputModule_.outputErrorLog();
// collect all data to I/O rank and assign to sol
Opm::data::Solution localCellData = {};
if (!isSubStep)
eclOutputModule_.assignToSolution(localCellData);
// add cell data to perforations for Rft output
if (!isSubStep)
eclOutputModule_.addRftDataToWells(localWellData, episodeIdx);
if (collectToIORank_.isParallel())
collectToIORank_.collect(localCellData, eclOutputModule_.getBlockData(), localWellData);
std::map<std::string, double> miscSummaryData;
std::map<std::string, std::vector<double>> regionData;
eclOutputModule_.outputFipLog(miscSummaryData, regionData, isSubStep);
// write output on I/O rank
if (collectToIORank_.isIORank()) {
const auto& eclState = simulator_.vanguard().eclState();
const auto& simConfig = eclState.getSimulationConfig();
// Add TCPU
if (totalCpuTime != 0.0)
miscSummaryData["TCPU"] = totalCpuTime;
bool enableDoublePrecisionOutput = EWOMS_GET_PARAM(TypeTag, bool, EclOutputDoublePrecision);
const Opm::data::Solution& cellData = collectToIORank_.isParallel() ? collectToIORank_.globalCellData() : localCellData;
const Opm::data::Wells& wellData = collectToIORank_.isParallel() ? collectToIORank_.globalWellData() : localWellData;
Opm::RestartValue restartValue(cellData, wellData);
const std::map<std::pair<std::string, int>, double>& blockData
= collectToIORank_.isParallel()
? collectToIORank_.globalBlockData()
: eclOutputModule_.getBlockData();
// Add suggested next timestep to extra data.
if (!isSubStep)
restartValue.addExtra("OPMEXTRA", std::vector<double>(1, nextStepSize));
if (simConfig.useThresholdPressure())
restartValue.addExtra("THRESHPR", Opm::UnitSystem::measure::pressure, simulator_.problem().thresholdPressure().data());
// first, create a tasklet to write the data for the current time step to disk
auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(*eclIO_,
episodeIdx,
isSubStep,
curTime,
restartValue,
miscSummaryData,
regionData,
blockData,
enableDoublePrecisionOutput);
// then, make sure that the previous I/O request has been completed and the
// number of incomplete tasklets does not increase between time steps
taskletRunner_->barrier();
// finally, start a new output writing job
taskletRunner_->dispatch(eclWriteTasklet);
}
}
void beginRestart()
{
bool enableHysteresis = simulator_.problem().materialLawManager()->enableHysteresis();
bool enableSwatinit = simulator_.vanguard().eclState().get3DProperties().hasDeckDoubleGridProperty("SWATINIT");
std::vector<Opm::RestartKey> solutionKeys{
{"PRESSURE", Opm::UnitSystem::measure::pressure},
{"SWAT", Opm::UnitSystem::measure::identity, static_cast<bool>(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))},
{"SGAS", Opm::UnitSystem::measure::identity, static_cast<bool>(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))},
{"TEMP" , Opm::UnitSystem::measure::temperature, enableEnergy},
{"RS", Opm::UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGas()},
{"RV", Opm::UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedOil()},
{"SOMAX", Opm::UnitSystem::measure::identity, simulator_.problem().vapparsActive()},
{"PCSWM_OW", Opm::UnitSystem::measure::identity, enableHysteresis},
{"KRNSW_OW", Opm::UnitSystem::measure::identity, enableHysteresis},
{"PCSWM_GO", Opm::UnitSystem::measure::identity, enableHysteresis},
{"KRNSW_GO", Opm::UnitSystem::measure::identity, enableHysteresis},
{"PPCW", Opm::UnitSystem::measure::pressure, enableSwatinit}
};
const auto& inputThpres = eclState().getSimulationConfig().getThresholdPressure();
std::vector<Opm::RestartKey> extraKeys = {{"OPMEXTRA", Opm::UnitSystem::measure::identity, false},
{"THRESHPR", Opm::UnitSystem::measure::pressure, inputThpres.active()}};
unsigned episodeIdx = simulator_.episodeIndex();
const auto& gridView = simulator_.vanguard().gridView();
unsigned numElements = gridView.size(/*codim=*/0);
eclOutputModule_.allocBuffers(numElements, episodeIdx, /*isSubStep=*/false, /*log=*/false);
auto restartValues = eclIO_->loadRestart(solutionKeys, extraKeys);
for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx);
eclOutputModule_.setRestart(restartValues.solution, elemIdx, globalIdx);
}
if (inputThpres.active()) {
Simulator& mutableSimulator = const_cast<Simulator&>(simulator_);
auto& thpres = mutableSimulator.problem().thresholdPressure();
const auto& thpresValues = restartValues.getExtra("THRESHPR");
thpres.setFromRestart(thpresValues);
}
restartTimeStepSize_ = restartValues.getExtra("OPMEXTRA")[0];
}
void endRestart()
{}
const EclOutputBlackOilModule<TypeTag>& eclOutputModule() const
{ return eclOutputModule_; }
Scalar restartTimeStepSize() const
{ return restartTimeStepSize_; }
private:
static bool enableEclOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
Opm::data::Solution computeTrans_() const
{
const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper();
const auto& cartDims = cartMapper.cartesianDimensions();
const int globalSize = cartDims[0]*cartDims[1]*cartDims[2];
Opm::data::CellData tranx = {Opm::UnitSystem::measure::transmissibility, std::vector<double>(globalSize), Opm::data::TargetType::INIT};
Opm::data::CellData trany = {Opm::UnitSystem::measure::transmissibility, std::vector<double>(globalSize), Opm::data::TargetType::INIT};
Opm::data::CellData tranz = {Opm::UnitSystem::measure::transmissibility, std::vector<double>(globalSize), Opm::data::TargetType::INIT};
for (size_t i = 0; i < tranx.data.size(); ++i) {
tranx.data[0] = 0.0;
trany.data[0] = 0.0;
tranz.data[0] = 0.0;
}
const auto& globalGridView = globalGrid_.leafGridView();
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> ElementMapper;
ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
#else
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout> ElementMapper;
ElementMapper globalElemMapper(globalGridView);
#endif
const auto& cartesianCellIdx = globalGrid_.globalCell();
const auto* globalTrans = &(simulator_.vanguard().globalTransmissibility());
if (!collectToIORank_.isParallel())
// in the sequential case we must use the transmissibilites defined by
// the problem. (because in the sequential case, the grid manager does
// not compute "global" transmissibilities for performance reasons. in
// the parallel case, the problem's transmissibilities can't be used
// because this object refers to the distributed grid and we need the
// sequential version here.)
globalTrans = &simulator_.problem().eclTransmissibilities();
auto elemIt = globalGridView.template begin</*codim=*/0>();
const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
const auto& elem = *elemIt;
auto isIt = globalGridView.ibegin(elem);
const auto& isEndIt = globalGridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
const auto& is = *isIt;
if (!is.neighbor())
continue; // intersection is on the domain boundary
unsigned c1 = globalElemMapper.index(is.inside());
unsigned c2 = globalElemMapper.index(is.outside());
if (c1 > c2)
continue; // we only need to handle each connection once, thank you.
int gc1 = std::min(cartesianCellIdx[c1], cartesianCellIdx[c2]);
int gc2 = std::max(cartesianCellIdx[c1], cartesianCellIdx[c2]);
if (gc2 - gc1 == 1)
tranx.data[gc1] = globalTrans->transmissibility(c1, c2);
if (gc2 - gc1 == cartDims[0])
trany.data[gc1] = globalTrans->transmissibility(c1, c2);
if (gc2 - gc1 == cartDims[0]*cartDims[1])
tranz.data[gc1] = globalTrans->transmissibility(c1, c2);
}
}
return {{"TRANX", tranx},
{"TRANY", trany},
{"TRANZ", tranz}};
}
Opm::NNC exportNncStructure_() const
{
Opm::NNC nnc = eclState().getInputNNC();
int nx = eclState().getInputGrid().getNX();
int ny = eclState().getInputGrid().getNY();
const auto& globalGridView = globalGrid_.leafGridView();
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> ElementMapper;
ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
#else
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout> ElementMapper;
ElementMapper globalElemMapper(globalGridView);
#endif
const auto* globalTrans = &(simulator_.vanguard().globalTransmissibility());
if (!collectToIORank_.isParallel()) {
// in the sequential case we must use the transmissibilites defined by
// the problem. (because in the sequential case, the grid manager does
// not compute "global" transmissibilities for performance reasons. in
// the parallel case, the problem's transmissibilities can't be used
// because this object refers to the distributed grid and we need the
// sequential version here.)
globalTrans = &simulator_.problem().eclTransmissibilities();
}
auto elemIt = globalGridView.template begin</*codim=*/0>();
const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
const auto& elem = *elemIt;
auto isIt = globalGridView.ibegin(elem);
const auto& isEndIt = globalGridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
const auto& is = *isIt;
if (!is.neighbor())
continue; // intersection is on the domain boundary
unsigned c1 = globalElemMapper.index(is.inside());
unsigned c2 = globalElemMapper.index(is.outside());
if (c1 > c2)
continue; // we only need to handle each connection once, thank you.
// TODO (?): use the cartesian index mapper to make this code work
// with grids other than Dune::CpGrid. The problem is that we need
// the a mapper for the sequential grid, not for the distributed one.
int cc1 = globalGrid_.globalCell()[c1];
int cc2 = globalGrid_.globalCell()[c2];
if (std::abs(cc1 - cc2) != 1 &&
std::abs(cc1 - cc2) != nx &&
std::abs(cc1 - cc2) != nx*ny)
{
nnc.addNNC(cc1, cc2, globalTrans->transmissibility(c1, c2));
}
}
}
return nnc;
}
struct EclWriteTasklet
: public TaskletInterface
{
Opm::EclipseIO& eclIO_;
int episodeIdx_;
bool isSubStep_;
double secondsElapsed_;
Opm::RestartValue restartValue_;
std::map<std::string, double> singleSummaryValues_;
std::map<std::string, std::vector<double>> regionSummaryValues_;
std::map<std::pair<std::string, int>, double> blockSummaryValues_;
bool writeDoublePrecision_;
explicit EclWriteTasklet(Opm::EclipseIO& eclIO,
int episodeIdx,
bool isSubStep,
double secondsElapsed,
Opm::RestartValue restartValue,
const std::map<std::string, double>& singleSummaryValues,
const std::map<std::string, std::vector<double>>& regionSummaryValues,
const std::map<std::pair<std::string, int>, double>& blockSummaryValues,
bool writeDoublePrecision)
: eclIO_(eclIO)
, episodeIdx_(episodeIdx)
, isSubStep_(isSubStep)
, secondsElapsed_(secondsElapsed)
, restartValue_(restartValue)
, singleSummaryValues_(singleSummaryValues)
, regionSummaryValues_(regionSummaryValues)
, blockSummaryValues_(blockSummaryValues)
, writeDoublePrecision_(writeDoublePrecision)
{ }
// callback to eclIO serial writeTimeStep method
void run()
{
eclIO_.writeTimeStep(episodeIdx_,
isSubStep_,
secondsElapsed_,
restartValue_,
singleSummaryValues_,
regionSummaryValues_,
blockSummaryValues_,
writeDoublePrecision_);
}
};
const Opm::EclipseState& eclState() const
{ return simulator_.vanguard().eclState(); }
const Simulator& simulator_;
CollectDataToIORankType collectToIORank_;
EclOutputBlackOilModule<TypeTag> eclOutputModule_;
std::unique_ptr<Opm::EclipseIO> eclIO_;
Grid globalGrid_;
std::unique_ptr<TaskletRunner> taskletRunner_;
Scalar restartTimeStepSize_;
};
} // namespace Ewoms
#endif

View File

@ -0,0 +1,999 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/**
* \file
*
* \brief Auxiliary routines that to solve the ODEs that emerge from the hydrostatic
* equilibrium problem
*/
#ifndef EWOMS_EQUILIBRATIONHELPERS_HH
#define EWOMS_EQUILIBRATIONHELPERS_HH
#include <opm/material/common/Tabulated1DFunction.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
// due to a bug in Equil.hpp, cstddef must be included before Equil.hpp
#include <cstddef>
#include <opm/parser/eclipse/EclipseState/InitConfig/Equil.hpp>
#include <memory>
/*
---- synopsis of EquilibrationHelpers.hpp ----
namespace Opm
{
namespace EQUIL {
namespace Miscibility {
class RsFunction;
class NoMixing;
template <class FluidSystem>
class RsVD;
template <class FluidSystem>
class RsSatAtContact;
}
class EquilReg;
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
struct PcEq;
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
double satFromPc(const MaterialLawManager& materialLawManager,
const int phase,
const int cell,
const double targetPc,
const bool increasing = false)
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
double satFromSumOfPcs(const MaterialLawManager& materialLawManager,
const int phase1,
const int phase2,
const int cell,
const double targetPc)
} // namespace Equil
} // namespace Opm
---- end of synopsis of EquilibrationHelpers.hpp ----
*/
namespace Ewoms {
/**
* Types and routines that collectively implement a basic
* ECLIPSE-style equilibration-based initialisation scheme.
*
* This namespace is intentionally nested to avoid name clashes
* with other parts of OPM.
*/
namespace EQUIL {
typedef Opm::BlackOilFluidSystem<double> FluidSystemSimple;
// Adjust oil pressure according to gas saturation and cap pressure
typedef Opm::SimpleModularFluidState<double,
/*numPhases=*/3,
/*numComponents=*/3,
FluidSystemSimple,
/*storePressure=*/false,
/*storeTemperature=*/false,
/*storeComposition=*/false,
/*storeFugacity=*/false,
/*storeSaturation=*/true,
/*storeDensity=*/false,
/*storeViscosity=*/false,
/*storeEnthalpy=*/false> SatOnlyFluidState;
/**
* Types and routines relating to phase mixing in
* equilibration calculations.
*/
namespace Miscibility {
/**
* Base class for phase mixing functions.
*/
class RsFunction
{
public:
virtual ~RsFunction() = default;
/**
* Function call operator.
*
* \param[in] depth Depth at which to calculate RS
* value.
*
* \param[in] press Pressure at which to calculate RS
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
virtual double operator()(const double depth,
const double press,
const double temp,
const double sat = 0.0) const = 0;
};
/**
* Type that implements "no phase mixing" policy.
*/
class NoMixing : public RsFunction
{
public:
virtual ~NoMixing() = default;
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RS
* value.
*
* \param[in] press Pressure at which to calculate RS
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. In "no mixing
* policy", this is identically zero.
*/
double
operator()(const double /* depth */,
const double /* press */,
const double /* temp */,
const double /* sat */ = 0.0) const
{
return 0.0;
}
};
/**
* Type that implements "dissolved gas-oil ratio"
* tabulated as a function of depth policy. Data
* typically taken from keyword 'RSVD'.
*/
template <class FluidSystem>
class RsVD : public RsFunction
{
public:
/**
* Constructor.
*
* \param[in] pvtRegionIdx The pvt region index
* \param[in] depth Depth nodes.
* \param[in] rs Dissolved gas-oil ratio at @c depth.
*/
RsVD(const int pvtRegionIdx,
const std::vector<double>& depth,
const std::vector<double>& rs)
: pvtRegionIdx_(pvtRegionIdx)
, rsVsDepth_(depth, rs)
{}
virtual ~RsVD() = default;
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RS
* value.
*
* \param[in] press Pressure at which to calculate RS
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
double operator()(const double depth,
const double press,
const double temp,
const double satGas = 0.0) const
{
if (satGas > 0.0) {
return satRs(press, temp);
}
else {
if (rsVsDepth_.xMin() > depth)
return rsVsDepth_.valueAt(0);
else if (rsVsDepth_.xMax() < depth)
return rsVsDepth_.valueAt(rsVsDepth_.numSamples() - 1);
return std::min(satRs(press, temp), rsVsDepth_.eval(depth, /*extrapolate=*/false));
}
}
private:
typedef Opm::Tabulated1DFunction<double> RsVsDepthFunc;
const int pvtRegionIdx_;
RsVsDepthFunc rsVsDepth_;
double satRs(const double press, const double temp) const
{
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
}
};
/**
* Type that implements "dissolved gas-oil ratio"
* tabulated as a function of depth policy. Data
* typically from keyword 'PBVD'.
*/
template <class FluidSystem>
class PBVD : public RsFunction
{
public:
/**
* Constructor.
*
* \param[in] pvtRegionIdx The pvt region index
* \param[in] depth Depth nodes.
* \param[in] pbub Bubble-point pressure at @c depth.
*/
PBVD(const int pvtRegionIdx,
const std::vector<double>& depth,
const std::vector<double>& pbub)
: pvtRegionIdx_(pvtRegionIdx)
, pbubVsDepth_(depth, pbub)
{}
virtual ~PBVD() = default;
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RS
* value.
*
* \param[in] Pressure in the cell
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
double operator()(const double depth,
const double cellPress,
const double temp,
const double satGas = 0.0) const
{
double press = cellPress;
if (satGas <= 0.0) {
if (pbubVsDepth_.xMin() > depth)
press = pbubVsDepth_.valueAt(0);
else if (pbubVsDepth_.xMax() < depth)
press = pbubVsDepth_.valueAt(pbubVsDepth_.numSamples() - 1);
else
press = pbubVsDepth_.eval(depth, /*extrapolate=*/false);
}
return satRs(std::min(press, cellPress), temp);
}
private:
typedef Opm::Tabulated1DFunction<double> PbubVsDepthFunc;
const int pvtRegionIdx_;
PbubVsDepthFunc pbubVsDepth_;
double satRs(const double press, const double temp) const
{
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
}
};
/**
* Type that implements "vaporized oil-gas ratio"
* tabulated as a function of depth policy. Data
* taken from keyword 'PDVD'.
*/
template <class FluidSystem>
class PDVD : public RsFunction
{
public:
/**
* Constructor.
*
* \param[in] pvtRegionIdx The pvt region index
* \param[in] depth Depth nodes.
* \param[in] pbub Dew-point pressure at @c depth.
*/
PDVD(const int pvtRegionIdx,
const std::vector<double>& depth,
const std::vector<double>& pdew)
: pvtRegionIdx_(pvtRegionIdx)
, pdewVsDepth_(depth, pdew)
{}
virtual ~PDVD() = default;
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RV
* value.
*
* \param[in] cellPress Pressure in the cell
*
* \param[in] temp Temperature at which to calculate RV
* value.
*
* \return Vaporized oil-gas ratio (RV) at depth @c
* depth and pressure @c press.
*/
double operator()(const double depth,
const double cellPress,
const double temp,
const double satOil = 0.0) const
{
double press = cellPress;
if (satOil <= 0.0) {
if (pdewVsDepth_.xMin() > depth)
press = pdewVsDepth_.valueAt(0);
else if (pdewVsDepth_.xMax() < depth)
press = pdewVsDepth_.valueAt(pdewVsDepth_.numSamples() - 1);
else
press = pdewVsDepth_.eval(depth, /*extrapolate=*/false);
}
return satRv(std::min(press, cellPress), temp);
}
private:
typedef Opm::Tabulated1DFunction<double> PdewVsDepthFunc;
const int pvtRegionIdx_;
PdewVsDepthFunc pdewVsDepth_;
double satRv(const double press, const double temp) const
{
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
}
};
/**
* Type that implements "vaporized oil-gas ratio"
* tabulated as a function of depth policy. Data
* typically taken from keyword 'RVVD'.
*/
template <class FluidSystem>
class RvVD : public RsFunction
{
public:
/**
* Constructor.
*
* \param[in] pvtRegionIdx The pvt region index
* \param[in] depth Depth nodes.
* \param[in] rv Dissolved gas-oil ratio at @c depth.
*/
RvVD(const int pvtRegionIdx,
const std::vector<double>& depth,
const std::vector<double>& rv)
: pvtRegionIdx_(pvtRegionIdx)
, rvVsDepth_(depth, rv)
{}
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RV
* value.
*
* \param[in] press Pressure at which to calculate RV
* value.
*
* \param[in] temp Temperature at which to calculate RV
* value.
*
* \return Vaporized oil-gas ratio (RV) at depth @c
* depth and pressure @c press.
*/
double operator()(const double depth,
const double press,
const double temp,
const double satOil = 0.0) const
{
if (std::abs(satOil) > 1e-16) {
return satRv(press, temp);
}
else {
if (rvVsDepth_.xMin() > depth)
return rvVsDepth_.valueAt(0);
else if (rvVsDepth_.xMax() < depth)
return rvVsDepth_.valueAt(rvVsDepth_.numSamples() - 1);
return std::min(satRv(press, temp), rvVsDepth_.eval(depth, /*extrapolate=*/false));
}
}
private:
typedef Opm::Tabulated1DFunction<double> RvVsDepthFunc;
const int pvtRegionIdx_;
RvVsDepthFunc rvVsDepth_;
double satRv(const double press, const double temp) const
{
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
}
};
/**
* Class that implements "dissolved gas-oil ratio" (Rs)
* as function of depth and pressure as follows:
*
* 1. The Rs at the gas-oil contact is equal to the
* saturated Rs value, RsSatContact.
*
* 2. The Rs elsewhere is equal to RsSatContact, but
* constrained to the saturated value as given by the
* local pressure.
*
* This should yield Rs-values that are constant below the
* contact, and decreasing above the contact.
*/
template <class FluidSystem>
class RsSatAtContact : public RsFunction
{
public:
/**
* Constructor.
*
* \param[in] pvtRegionIdx The pvt region index
* \param[in] pContact oil pressure at the contact
* \param[in] T_contact temperature at the contact
*/
RsSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact)
: pvtRegionIdx_(pvtRegionIdx)
{
rsSatContact_ = satRs(pContact, T_contact);
}
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RS
* value.
*
* \param[in] press Pressure at which to calculate RS
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
double operator()(const double /* depth */,
const double press,
const double temp,
const double satGas = 0.0) const
{
if (satGas > 0.0) {
return satRs(press, temp);
}
else {
return std::min(satRs(press, temp), rsSatContact_);
}
}
private:
const int pvtRegionIdx_;
double rsSatContact_;
double satRs(const double press, const double temp) const
{
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
}
};
/**
* Class that implements "vaporized oil-gas ratio" (Rv)
* as function of depth and pressure as follows:
*
* 1. The Rv at the gas-oil contact is equal to the
* saturated Rv value, RvSatContact.
*
* 2. The Rv elsewhere is equal to RvSatContact, but
* constrained to the saturated value as given by the
* local pressure.
*
* This should yield Rv-values that are constant below the
* contact, and decreasing above the contact.
*/
template <class FluidSystem>
class RvSatAtContact : public RsFunction
{
public:
/**
* Constructor.
*
* \param[in] pvtRegionIdx The pvt region index
* \param[in] pContact oil pressure at the contact
* \param[in] T_contact temperature at the contact
*/
RvSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact)
:pvtRegionIdx_(pvtRegionIdx)
{
rvSatContact_ = satRv(pContact, T_contact);
}
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RV
* value.
*
* \param[in] press Pressure at which to calculate RV
* value.
*
* \param[in] temp Temperature at which to calculate RV
* value.
*
* \return Dissolved oil-gas ratio (RV) at depth @c
* depth and pressure @c press.
*/
double operator()(const double /*depth*/,
const double press,
const double temp,
const double satOil = 0.0) const
{
if (satOil > 0.0) {
return satRv(press, temp);
}
else {
return std::min(satRv(press, temp), rvSatContact_);
}
}
private:
const int pvtRegionIdx_;
double rvSatContact_;
double satRv(const double press, const double temp) const
{
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
}
};
} // namespace Miscibility
/**
* Aggregate information base of an equilibration region.
*
* Provides inquiry methods for retrieving depths of contacs
* and pressure values as well as a means of calculating fluid
* densities, dissolved gas-oil ratio and vapourised oil-gas
* ratios.
*
* \tparam DensCalc Type that provides access to a phase
* density calculation facility. Must implement an operator()
* declared as
* <CODE>
* std::vector<double>
* operator()(const double press,
* const std::vector<double>& svol)
* </CODE>
* that calculates the phase densities of all phases in @c
* svol at fluid pressure @c press.
*/
class EquilReg
{
public:
/**
* Constructor.
*
* \param[in] rec Equilibration data of current region.
* \param[in] rs Calculator of dissolved gas-oil ratio.
* \param[in] rv Calculator of vapourised oil-gas ratio.
* \param[in] pvtRegionIdx The pvt region index
*/
EquilReg(const Opm::EquilRecord& rec,
std::shared_ptr<Miscibility::RsFunction> rs,
std::shared_ptr<Miscibility::RsFunction> rv,
const int pvtIdx)
: rec_ (rec)
, rs_ (rs)
, rv_ (rv)
, pvtIdx_ (pvtIdx)
{}
/**
* Type of dissolved gas-oil ratio calculator.
*/
typedef Miscibility::RsFunction CalcDissolution;
/**
* Type of vapourised oil-gas ratio calculator.
*/
typedef Miscibility::RsFunction CalcEvaporation;
/**
* Datum depth in current region
*/
double datum() const { return this->rec_.datumDepth(); }
/**
* Pressure at datum depth in current region.
*/
double pressure() const { return this->rec_.datumDepthPressure(); }
/**
* Depth of water-oil contact.
*/
double zwoc() const { return this->rec_.waterOilContactDepth(); }
/**
* water-oil capillary pressure at water-oil contact.
*
* \return P_o - P_w at WOC.
*/
double pcowWoc() const { return this->rec_.waterOilContactCapillaryPressure(); }
/**
* Depth of gas-oil contact.
*/
double zgoc() const { return this->rec_.gasOilContactDepth(); }
/**
* Gas-oil capillary pressure at gas-oil contact.
*
* \return P_g - P_o at GOC.
*/
double pcgoGoc() const { return this->rec_.gasOilContactCapillaryPressure(); }
/**
* Retrieve dissolved gas-oil ratio calculator of current
* region.
*/
const CalcDissolution&
dissolutionCalculator() const { return *this->rs_; }
/**
* Retrieve vapourised oil-gas ratio calculator of current
* region.
*/
const CalcEvaporation&
evaporationCalculator() const { return *this->rv_; }
/**
* Retrieve pvtIdx of the region.
*/
int pvtIdx() const { return this->pvtIdx_; }
private:
Opm::EquilRecord rec_; /**< Equilibration data */
std::shared_ptr<Miscibility::RsFunction> rs_; /**< RS calculator */
std::shared_ptr<Miscibility::RsFunction> rv_; /**< RV calculator */
const int pvtIdx_;
};
/// Functor for inverting capillary pressure function.
/// Function represented is
/// f(s) = pc(s) - targetPc
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
struct PcEq
{
PcEq(const MaterialLawManager& materialLawManager,
const int phase,
const int cell,
const double targetPc)
: materialLawManager_(materialLawManager),
phase_(phase),
cell_(cell),
targetPc_(targetPc)
{}
double operator()(double s) const
{
const auto& matParams = materialLawManager_.materialLawParams(cell_);
SatOnlyFluidState fluidState;
fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
fluidState.setSaturation(phase_, s);
double pc[FluidSystem::numPhases];
std::fill(pc, pc + FluidSystem::numPhases, 0.0);
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_];
return pcPhase - targetPc_;
}
private:
const MaterialLawManager& materialLawManager_;
const int phase_;
const int cell_;
const double targetPc_;
};
template <class FluidSystem, class MaterialLawManager>
double minSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell)
{
const auto& scaledDrainageInfo =
materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
// Find minimum and maximum saturations.
switch(phase) {
case FluidSystem::waterPhaseIdx:
return scaledDrainageInfo.Swl;
case FluidSystem::gasPhaseIdx:
return scaledDrainageInfo.Sgl;
case FluidSystem::oilPhaseIdx:
throw std::runtime_error("Min saturation not implemented for oil phase.");
default:
throw std::runtime_error("Unknown phaseIdx .");
}
return -1.0;
}
template <class FluidSystem, class MaterialLawManager>
double maxSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell)
{
const auto& scaledDrainageInfo =
materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
// Find minimum and maximum saturations.
switch(phase) {
case FluidSystem::waterPhaseIdx:
return scaledDrainageInfo.Swu;
case FluidSystem::gasPhaseIdx:
return scaledDrainageInfo.Sgu;
case FluidSystem::oilPhaseIdx:
throw std::runtime_error("Max saturation not implemented for oil phase.");
default:
throw std::runtime_error("Unknown phaseIdx .");
}
return -1.0;
}
/// Compute saturation of some phase corresponding to a given
/// capillary pressure.
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
double satFromPc(const MaterialLawManager& materialLawManager,
const int phase,
const int cell,
const double targetPc,
const bool increasing = false)
{
// Find minimum and maximum saturations.
double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
// Create the equation f(s) = pc(s) - targetPc
const PcEq<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase, cell, targetPc);
double f0 = f(s0);
double f1 = f(s1);
if (f0 <= 0.0)
return s0;
else if (f1 >= 0.0)
return s1;
assert(f0 > 0 && f1 < 0);
const int maxIter = 60;
const double tol = 1e-10;
// regula falsi with the "Pegasus" method to avoid stagnation
for (int iter = 0; iter < maxIter; ++ iter) {
// determine the pivot
double s = (s1*f0 - s0*f1)/(f0 - f1);
// adapt the interval
double fs = f(s);
if (fs == 0.0)
return s;
else if ((fs > 0.0) == (f0 > 0.0)) {
// update interval and reverse
s0 = s1;
f0 = f1;
}
else
// "Pegasus" method
f0 *= f1/(f1 + fs);
s1 = s;
f1 = fs;
// check for convergence
if (std::abs(s1 - s0) < tol)
return (s1 + s0)/2;
}
throw std::runtime_error("Could not find solution for PcEq = 0.0 after "+std::to_string(maxIter)
+" iterations.");
}
/// Functor for inverting a sum of capillary pressure functions.
/// Function represented is
/// f(s) = pc1(s) + pc2(1 - s) - targetPc
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
struct PcEqSum
{
PcEqSum(const MaterialLawManager& materialLawManager,
const int phase1,
const int phase2,
const int cell,
const double targetPc)
: materialLawManager_(materialLawManager),
phase1_(phase1),
phase2_(phase2),
cell_(cell),
targetPc_(targetPc)
{}
double operator()(double s) const
{
const auto& matParams = materialLawManager_.materialLawParams(cell_);
SatOnlyFluidState fluidState;
fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
fluidState.setSaturation(phase1_, s);
fluidState.setSaturation(phase2_, 1.0 - s);
double pc[FluidSystem::numPhases];
std::fill(pc, pc + FluidSystem::numPhases, 0.0);
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
double sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
double pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_];
double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
double pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_];
return pc1 + pc2 - targetPc_;
}
private:
const MaterialLawManager& materialLawManager_;
const int phase1_;
const int phase2_;
const int cell_;
const double targetPc_;
};
/// Compute saturation of some phase corresponding to a given
/// capillary pressure, where the capillary pressure function
/// is given as a sum of two other functions.
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
double satFromSumOfPcs(const MaterialLawManager& materialLawManager,
const int phase1,
const int phase2,
const int cell,
const double targetPc)
{
// Find minimum and maximum saturations.
double s0 = minSaturations<FluidSystem>(materialLawManager, phase1, cell);
double s1 = maxSaturations<FluidSystem>(materialLawManager, phase1, cell);
// Create the equation f(s) = pc1(s) + pc2(1-s) - targetPc
const PcEqSum<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase1, phase2, cell, targetPc);
double f0 = f(s0);
double f1 = f(s1);
if (f0 <= 0.0)
return s0;
else if (f1 >= 0.0)
return s1;
assert(f0 > 0.0 && f1 < 0.0);
const int maxIter = 60;
const double tol = 1e-10;
// regula falsi with the "Pegasus" method to avoid stagnation
for (int iter = 0; iter < maxIter; ++ iter) {
// determine the pivot
double s = (s1*f0 - s0*f1)/(f0 - f1);
// adapt the interval
double fs = f(s);
if (fs == 0.0)
return s;
else if ((fs > 0.0) == (f0 > 0.0)) {
// update interval and reverse
s0 = s1;
f0 = f1;
}
else
// "Pegasus" method
f0 *= f1 / (f1 + fs);
s1 = s;
f1 = fs;
// check for convergence
if (std::abs(s1 - s0) < tol)
return (s1 + s0)/2;
}
throw std::runtime_error("Could not find solution for PcEqSum = 0.0 after "+std::to_string(maxIter)
+" iterations.");
}
/// Compute saturation from depth. Used for constant capillary pressure function
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
double satFromDepth(const MaterialLawManager& materialLawManager,
const double cellDepth,
const double contactDepth,
const int phase,
const int cell,
const bool increasing = false)
{
const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
if (cellDepth < contactDepth) {
return s0;
}
else {
return s1;
}
}
/// Return true if capillary pressure function is constant
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
bool isConstPc(const MaterialLawManager& materialLawManager,
const int phase,
const int cell)
{
// Create the equation f(s) = pc(s);
const PcEq<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase, cell, 0);
const double f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell));
const double f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell));
return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
}
} // namespace Equil
} // namespace Ewoms
#endif // EWOMS_EQUILIBRATIONHELPERS_HH

1188
ebos/equil/initstateequil.hh Normal file

File diff suppressed because it is too large Load Diff

217
ebos/equil/regionmapping.hh Normal file
View File

@ -0,0 +1,217 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef OPM_REGIONMAPPING_HH
#define OPM_REGIONMAPPING_HH
#include <unordered_map>
#include <vector>
namespace Ewoms {
/**
* Forward and reverse mappings between cells and
* regions/partitions (e.g., the ECLIPSE-style 'SATNUM',
* 'PVTNUM', or 'EQUILNUM' arrays).
*
* \tparam Region Type of a forward region mapping. Expected
* to provide indexed access through
* operator[]() as well as inner types
* 'valueType', 'size_type', and
* 'const_iterator'.
*/
template <class Region = std::vector<int>>
class RegionMapping
{
public:
/**
* Constructor.
*
* \param[in] reg Forward region mapping, restricted to
* active cells only.
*/
explicit
RegionMapping(const Region& reg)
: reg_(reg)
{
rev_.init(reg_);
}
/**
* Type of forward (cell-to-region) mapping result.
* Expected to be an integer.
*/
typedef typename Region::value_type RegionId;
/**
* Type of reverse (region-to-cell) mapping (element)
* result.
*/
typedef typename Region::size_type CellId;
/**
* Type of reverse region-to-cell range bounds and
* iterators.
*/
typedef typename std::vector<CellId>::const_iterator CellIter;
class Range
{
public:
typedef CellIter iterator;
typedef CellIter const_iterator;
Range()
{};
Range(const CellIter& beg, const CellIter& en)
: begin_(beg)
, end_(en)
{};
Range(const Range&) = default;
CellIter& begin() { return begin_; }
const CellIter& begin() const { return begin_; }
const CellIter& end() const { return end_; }
bool empty() const
{ return begin_ == end_; }
size_t size() const
{
size_t ret = 0;
for (CellIter it = begin(); it != end(); ++it)
++ret;
return ret;
}
private:
CellIter begin_;
CellIter end_;
};
/**
* Compute region number of given active cell.
*
* \param[in] c Active cell
* \return Region to which @c c belongs.
*/
RegionId region(const CellId c) const
{ return reg_[c]; }
const std::vector<RegionId>&
activeRegions() const
{
return rev_.active;
}
/**
* Extract active cells in particular region.
*
* \param[in] r Region number
*
* \return Range of active cells in region @c r. Empty if @c r is
* not an active region.
*/
Range cells(const RegionId r) const
{
const auto id = rev_.binid.find(r);
if (id == rev_.binid.end()) {
// Region 'r' not an active region. Return empty.
return Range(rev_.c.end(), rev_.c.end());
}
const auto i = id->second;
return Range(rev_.c.begin() + rev_.p[i + 0],
rev_.c.begin() + rev_.p[i + 1]);
}
private:
/**
* Copy of forward region mapping (cell-to-region).
*/
Region reg_;
/**
* Reverse mapping (region-to-cell).
*/
struct {
typedef typename std::vector<CellId>::size_type Pos;
std::unordered_map<RegionId, Pos> binid;
std::vector<RegionId> active;
std::vector<Pos> p; /**< Region start pointers */
std::vector<CellId> c; /**< Region cells */
/**
* Compute reverse mapping. Standard linear insertion
* sort algorithm.
*/
void
init(const Region& reg)
{
binid.clear();
for (const auto& r : reg) {
++binid[r];
}
p .clear(); p.emplace_back(0);
active.clear();
{
Pos n = 0;
for (auto& id : binid) {
active.push_back(id.first);
p .push_back(id.second);
id.second = n++;
}
}
for (decltype(p.size()) i = 1, sz = p.size(); i < sz; ++i) {
p[0] += p[i];
p[i] = p[0] - p[i];
}
assert (p[0] == static_cast<Pos>(reg.size()));
c.resize(reg.size());
{
CellId i = 0;
for (const auto& r : reg) {
auto& pos = p[binid[r] + 1];
c[pos++] = i++;
}
}
p[0] = 0;
}
} rev_; /**< Reverse mapping instance */
};
} // namespace Ewoms
#endif // OPM_REGIONMAPPING_HH

93
ebos/femcpgridcompat.hh Normal file
View File

@ -0,0 +1,93 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \brief This file ensures that ebos can be compiled in the presence of dune-fem
*
* It implements a few work-arounds for some incompatibilities with the Dune grid
* interface of CpGrid. A better way would be to make CpGrid conforming.
*/
#ifndef EWOMS_FEM_CPGRID_COMPAT_HH
#define EWOMS_FEM_CPGRID_COMPAT_HH
#if HAVE_DUNE_FEM
#include <dune/common/version.hh>
#include <dune/fem/gridpart/common/gridpart.hh>
#include <dune/fem/misc/compatibility.hh>
#include <dune/fem/io/streams/streams.hh>
namespace Dune {
namespace cpgrid {
template <int codim>
class Entity;
template <int codim>
class EntityPointer;
}
#if DUNE_VERSION_NEWER(DUNE_FEM, 2, 6)
template <int dim, int cdim>
auto referenceElement(const Dune::cpgrid::Geometry<dim, cdim>& geo)
-> decltype(referenceElement<double, dim>(geo.type()))
{ return referenceElement<double, dim>(geo.type()); }
#endif
// specialization of dune-fem compatiblity functions for CpGrid, since CpGrid does not use the interface classes.
namespace Fem {
////////////////////////////////////////////////////////////
//
// make_entity for CpGrid entities
//
////////////////////////////////////////////////////////////
template <int codim>
inline Dune::cpgrid::Entity<codim> make_entity(const Dune::cpgrid::EntityPointer<codim>& entityPointer)
{ return *entityPointer; }
template <int codim>
inline Dune::cpgrid::Entity<codim> make_entity(Dune::cpgrid::Entity<codim> entity)
{ return std::move(entity); }
////////////////////////////////////////////////////////////
//
// GridEntityAccess for CpGrid entities
//
////////////////////////////////////////////////////////////
template<int codim>
struct GridEntityAccess<Dune::cpgrid::Entity<codim> >
{
typedef Dune::cpgrid::Entity<codim> EntityType;
typedef EntityType GridEntityType;
static const GridEntityType& gridEntity(const EntityType& entity)
{ return entity; }
};
} // namespace Fem
} // end namespace Dune
#endif // #if HAVE_DUNE_FEM
#endif // EWOMS_FEM_CPGRID_COMPAT_HH

65
ebos/tracervdtable.hh Normal file
View File

@ -0,0 +1,65 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Ewoms::TracerVdTable
*/
#ifndef EWOMS_TRACER_VD_TABLE_HH
#define EWOMS_TRACER_VD_TABLE_HH
#include <opm/parser/eclipse/EclipseState/Tables/SimpleTable.hpp>
namespace Ewoms {
/*!
* \brief A class that contains tracer concentration vs depth table
*/
class TracerVdTable : public Opm::SimpleTable
{
public:
TracerVdTable(const Opm::DeckItem& item)
{
this->m_schema.addColumn(Opm::ColumnSchema("DEPTH", Opm::Table::STRICTLY_INCREASING, Opm::Table::DEFAULT_NONE));
this->m_schema.addColumn(Opm::ColumnSchema("TRACER_CONCENTRATION", Opm::Table::RANDOM, Opm::Table::DEFAULT_NONE));
Opm::SimpleTable::init(item);
}
/*!
* \brief Return the depth column
*/
const Opm::TableColumn& getDepthColumn() const
{ return Opm::SimpleTable::getColumn(0); }
/*!
* \brief Return the tracer concentration column
*/
const Opm::TableColumn& getTracerConcentration() const
{ return Opm::SimpleTable::getColumn(1); }
};
}
#endif // TRACERVDTABLE_HH

170
ebos/vtkecltracermodule.hh Normal file
View File

@ -0,0 +1,170 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
* \copydoc Ewoms::VtkEclTracerModule
*/
#ifndef EWOMS_VTK_ECL_TRACER_MODULE_HH
#define EWOMS_VTK_ECL_TRACER_MODULE_HH
#include <ewoms/io/vtkmultiwriter.hh>
#include <ewoms/io/baseoutputmodule.hh>
#include <ewoms/common/propertysystem.hh>
#include <ewoms/common/parametersystem.hh>
#include <ewoms/models/blackoil/blackoilproperties.hh>
#include <dune/common/fvector.hh>
#include <cstdio>
BEGIN_PROPERTIES
// create new type tag for the VTK tracer output
NEW_TYPE_TAG(VtkEclTracer);
// create the property tags needed for the tracer model
NEW_PROP_TAG(EnableVtkOutput);
NEW_PROP_TAG(VtkOutputFormat);
NEW_PROP_TAG(VtkWriteEclTracerConcentration);
// set default values for what quantities to output
SET_BOOL_PROP(VtkEclTracer, VtkWriteEclTracerConcentration, false);
END_PROPERTIES
namespace Ewoms {
/*!
* \ingroup Vtk
*
* \brief VTK output module for the tracer model's parameters.
*/
template <class TypeTag>
class VtkEclTracerModule : public BaseOutputModule<TypeTag>
{
typedef BaseOutputModule<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat);
typedef Ewoms::VtkMultiWriter<GridView, vtkFormat> VtkMultiWriter;
typedef typename ParentType::ScalarBuffer ScalarBuffer;
public:
VtkEclTracerModule(const Simulator& simulator)
: ParentType(simulator)
{ }
/*!
* \brief Register all run-time parameters for the tracer VTK output
* module.
*/
static void registerParameters()
{
EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteEclTracerConcentration,
"Include the tracer concentration "
"in the VTK output files");
}
/*!
* \brief Allocate memory for the scalar fields we would like to
* write to the VTK file.
*/
void allocBuffers()
{
if (eclTracerConcentrationOutput_()){
const auto& tracerModel = this->simulator_.problem().tracerModel();
eclTracerConcentration_.resize(tracerModel.numTracers());
for(size_t tracerIdx=0; tracerIdx<eclTracerConcentration_.size();++tracerIdx){
this->resizeScalarBuffer_(eclTracerConcentration_[tracerIdx]);
}
}
}
/*!
* \brief Modify the internal buffers according to the intensive quantities relevant for
* an element
*/
void processElement(const ElementContext& elemCtx)
{
if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput))
return;
const auto& tracerModel = elemCtx.problem().tracerModel();
for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
if (eclTracerConcentrationOutput_()){
for(size_t tracerIdx=0; tracerIdx<eclTracerConcentration_.size();++tracerIdx){
eclTracerConcentration_[tracerIdx][globalDofIdx] = tracerModel.tracerConcentration(tracerIdx, globalDofIdx);
}
}
}
}
/*!
* \brief Add all buffers to the VTK output writer.
*/
void commitBuffers(BaseOutputWriter& baseWriter)
{
VtkMultiWriter *vtkWriter = dynamic_cast<VtkMultiWriter*>(&baseWriter);
if (!vtkWriter)
return;
if (eclTracerConcentrationOutput_()){
const auto& tracerModel = this->simulator_.problem().tracerModel();
for(size_t tracerIdx=0; tracerIdx<eclTracerConcentration_.size();++tracerIdx){
const std::string tmp = "tracerConcentration_" + tracerModel.tracerName(tracerIdx);
this->commitScalarBuffer_(baseWriter,tmp.c_str(), eclTracerConcentration_[tracerIdx]);
}
}
}
private:
static bool eclTracerConcentrationOutput_()
{
static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteEclTracerConcentration);
return val;
}
std::vector<ScalarBuffer> eclTracerConcentration_;
};
} // namespace Ewoms
#endif

View File

@ -0,0 +1,658 @@
-- Synthetic test deck based on Norne. This data set is meant to be a simple,
-- well-documented deck for the behaviour of SUMMARY specified output. Data
-- is mostly entered to *traceable* and does not necessarily make sense from
-- a simulation point of view.
START
10 MAI 2007 /
RUNSPEC
TITLE
SUMMARYTESTS
-- A simple 10x10x10 cube. Simple to reason about, large enough for all tests
DIMENS
10 10 10 /
REGDIMS
3 /
OIL
GAS
WATER
DISGAS
METRIC
GRID
DX
1000*1 /
DY
1000*1 /
DZ
1000*1 /
TOPS
100*1 /
ACTNUM
1000*1/
PORO
500*0.1 500*0.2/
PERMX
1000*500 /
REGIONS
FIPNUM
400*1
200*2
400*3 /
PROPS
PVTW
1 1000 0 0.318 0.0 /
ROCK
14.7 3E-6 /
SWOF
0.12 0 1 0
0.18 4.64876033057851E-008 1 0
0.24 0.000000186 0.997 0
0.3 4.18388429752066E-007 0.98 0
0.36 7.43801652892562E-007 0.7 0
0.42 1.16219008264463E-006 0.35 0
0.48 1.67355371900826E-006 0.2 0
0.54 2.27789256198347E-006 0.09 0
0.6 2.97520661157025E-006 0.021 0
0.66 3.7654958677686E-006 0.01 0
0.72 4.64876033057851E-006 0.001 0
0.78 0.000005625 0.0001 0
0.84 6.69421487603306E-006 0 0
0.91 8.05914256198347E-006 0 0
1 0.00001 0 0 /
SGOF
0 0 1 0
0.001 0 1 0
0.02 0 0.997 0
0.05 0.005 0.980 0
0.12 0.025 0.700 0
0.2 0.075 0.350 0
0.25 0.125 0.200 0
0.3 0.190 0.090 0
0.4 0.410 0.021 0
0.45 0.60 0.010 0
0.5 0.72 0.001 0
0.6 0.87 0.0001 0
0.7 0.94 0.000 0
0.85 0.98 0.000 0
0.88 0.984 0.000 0 /
DENSITY
53.66 64.49 0.0533 /
PVDG
1 100 1
10 10 1 /
PVTO
0 1 10 1 /
1 1 10.001 1
10 1 1 /
/
SUMMARY
DATE
PERFORMA
--
-- Field Data
-- Production Rates
FVPR
FWPR
FWPRH
FOPR
FOPRH
FGPR
FGPRH
FLPR
FLPRH
FGSR
FGCR
FNPR -- solvent
--FTPRSEA
-- Injection Rates
FVIR
FWIR
FWIRH
FGIR
FNIR -- solvent
FGIRH
-- Production Cummulatives
FVPT
FWPT
FOPT
FLPT
FLPTH
FGPT
FNPT
FOPTH
FGPTH
FWPTH
FGST
FGCT
-- Injection Cummulatives
FVIT
FWIT
FWITH
FGIT
FNIT
FGITH
-- In place
FWIP
FOIP
FGIP
-- Ratios
FWCT
FWCTH
FGOR
FGORH
-- From model2
FMWPR
FMWIN
FOE
-- Pressures
FPR
FPRP
BPR
1 1 1 /
1 1 2 /
1 1 3 /
1 1 4 /
1 1 5 /
1 1 6 /
1 1 7 /
1 1 8 /
1 1 9 /
1 1 10 /
2 1 10 / -- This cell is not ACTIVE
/
BSGAS
1 1 1 /
/
BSWAT
1 1 1 /
/
-- Region data
RPR
/
ROPT
/
RGPT
/
RWPT
/
RGFT
/
RWFT
/
ROIP
/
ROP
/
ROPR
/
RGPR
/
RWPR
/
RGIR
/
RGIT
/
RWIR
/
RWIT
/
RWPT
/
ROIPL
/
ROIPG
/
RGIP
/
RGIPL
/
RGIPG
/
RWIP
/
RPPO
/
-- Group data --
GPR
/
GLPR
/
GOPT
/
GGPT
/
GWPT
/
GNPT
/
GOPR
/
GGPR
/
GWPR
/
GWPRH
/
GGIR
/
GNPR
/
GNIR
/
GGIRH
/
GGIT
/
GNIT
/
GGITH
/
GWCT
/
GWCTH
/
GGOR
/
GGORH
/
GWIR
/
GWIT
/
GWIRH
/
GWITH
/
GOPRH
/
GGPRH
/
GLPRH
/
GWPTH
/
GOPTH
/
GGPTH
/
GLPTH
/
GPRG
/
GPRW
/
GOPTF
/
GOPTS
/
GOPTH
/
GOPRF
/
GOPRS
/
GOPRH
/
GGPTF
/
GGPTS
/
GGPTH
/
GGPTF
/
GGPTS
/
GGPTH
/
GGLR
/
GGLIR
/
GGLRH
/
GVPR
/
GVPT
/
GMCTP
/
GOPP
/
GVIR
/
GVIT
/
GVPRT
/
GMWPR
/
GMWIN
/
-- Well Data
-- Production Rates
WWPR
/
WWPRH
/
WOPR
/
WOPRH
/
WGPR
/
WNPR
/
WGPRH
/
WLPR
/
WLPRH
/
WLPT
/
WLPTH
/
-- Injection Rates
WWIR
W_3
/
WWIT
W_3
/
WWIRH
W_3
/
WWITH
W_3
/
WGIT
W_3
/
WGIR
W_3
/
WGIRH
W_3
/
WGITH
W_3
/
WNIR
W_3
/
WNIT
W_3
/
-- Production Cummulatives
WWPT
/
WWPTH
/
WOPT
/
WOPTH
/
WGPT
/
WGPTH
/
WNPT
/
-- Tracers
--WTPRSEA
--/
--WTPTSEA
--/
-- Injection Cummulatives
WWIT
W_3
/
-- Ratios
WWCT
/
WWCTH
/
WGOR
/
WGORH
/
WGLR
/
WGLRH
/
-- Performance
WBHP
/
WBHPH
/
WTHP
/
WTHPH
/
WPI
/
WBP
/
WBP4
/
-- from model2
WOPTF
/
WOPTS
/
WOPTH
/
WOPRS
/
WOPRF
/
WGPTF
/
WGPTS
/
WGPRF
/
WTPRS
/
WGLIR
/
WVPR
/
WVPT
/
WOPP
/
WVIR
/
WVIT
/
WMCTL
/
-- Water injection per connection
CWIR
* /
/
-- Gas injection on 3 1 1 (45)
CGIR
'W_3' 3 1 1 /
/
CWIT
'W_3' /
/
CGIT
* /
/
-- Production per connection
-- Using all the different ways of specifying connections here
-- as an informal test that we still get the data we want
CWPR
'W_1' 1 1 1 /
/
COPR
'W_1' /
'W_2' /
'W_3' /
/
CGPR
'*' /
/
CNFR
'*' /
/
CNPT
'*' /
/
CNIT
'*' /
/
CWPT
'W_1' 1 1 1 /
/
COPT
'W_1' /
/
CGPT
'W_1' /
'W_2' /
'W_3' /
/
---- Connection production rates
----CGFR
----'E-4AH' /
----/
----CWFR
----'E-2H' /
----/
SOLUTION
SWAT
1000*0.2 /
SGAS
1000*0.0 /
PRESSURE
100*1.0
100*2.0
100*3.0
100*4.0
100*5.0
100*6.0
100*7.0
100*8.0
100*9.0
100*10.0/
RS
1000*0 /
SCHEDULE
-- Three wells, two producers (so that we can form a group) and one injector
WELSPECS
'W_1' 'G_1' 1 1 3.33 'OIL' 7* /
'W_2' 'G_1' 2 1 3.33 'OIL' 7* /
'W_3' 'G_2' 3 1 3.92 'WATER' 7* /
/
-- Completion data.
COMPDAT
-- 'Well' I J K1 K2
-- Passing 0 to I/J means they'll get the well head I/J
W_1 0 0 1 1 / -- Active index: 0
W_2 0 0 1 1 / -- Active index: 1
W_2 0 0 2 2 / -- Active index: 101
W_3 0 0 1 1 / -- Active index: 2
/
WCONHIST
-- history rates are set so that W_1 produces 1, W_2 produces 2 etc.
-- index.offset.
-- organised as oil-water-gas
W_1 SHUT ORAT 10.1 10 10.2 /
W_2 SHUT ORAT 20.1 20 20.2 /
/
WCONINJH
-- Injection historical rates (water only, as we only support pure injectors)
W_3 WATER STOP 30.0 /
/
TSTEP
-- register time steps (in days). This allows us to write *two* report steps (1
-- and 2. Without this, totals/accumulations would fail (segfault) when looking
-- up historical rates and volumes. These volumes however don't change, i.e.
-- every time step has the same set of values
10 10 /
-- Register a fourth well with completions later. This ensure we handle when
-- wells are registered or activated later in a simulation
WELSPECS
'W_4' 'G_3' 1 1 3.33 'OIL' 7* /
/
COMPDAT
W_4 1 1 3 3 /
/
TSTEP
10 10 /

115
tests/equil_base.DATA Normal file
View File

@ -0,0 +1,115 @@
RUNSPEC
WATER
GAS
OIL
METRIC
DIMENS
10 1 10 /
GRID
DX
100*1 /
DY
100*1 /
DZ
100*1 /
TOPS
10*0. /
PORO
100*0.3 /
PERMX
100*500 /
PROPS
PVTW
4017.55 1.038 3.22E-6 0.318 0.0 /
ROCK
14.7 3E-6 /
SWOF
0.12 0 1 0
0.18 4.64876033057851E-008 1 0
0.24 0.000000186 0.997 0
0.3 4.18388429752066E-007 0.98 0
0.36 7.43801652892562E-007 0.7 0
0.42 1.16219008264463E-006 0.35 0
0.48 1.67355371900826E-006 0.2 0
0.54 2.27789256198347E-006 0.09 0
0.6 2.97520661157025E-006 0.021 0
0.66 3.7654958677686E-006 0.01 0
0.72 4.64876033057851E-006 0.001 0
0.78 0.000005625 0.0001 0
0.84 6.69421487603306E-006 0 0
0.91 8.05914256198347E-006 0 0
1 0.00001 0 0 /
SGOF
0 0 1 0
0.001 0 1 0
0.02 0 0.997 0
0.05 0.005 0.980 0
0.12 0.025 0.700 0
0.2 0.075 0.350 0
0.25 0.125 0.200 0
0.3 0.190 0.090 0
0.4 0.410 0.021 0
0.45 0.60 0.010 0
0.5 0.72 0.001 0
0.6 0.87 0.0001 0
0.7 0.94 0.000 0
0.85 0.98 0.000 0
0.88 0.984 0.000 0 /
DENSITY
53.66 64.49 0.0533 /
PVDG
14.700 166.666 0.008000
264.70 12.0930 0.009600
514.70 6.27400 0.011200
1014.7 3.19700 0.014000
2014.7 1.61400 0.018900
2514.7 1.29400 0.020800
3014.7 1.08000 0.022800
4014.7 0.81100 0.026800
5014.7 0.64900 0.030900
9014.7 0.38600 0.047000 /
PVTO
0.0010 14.7 1.0620 1.0400 /
0.0905 264.7 1.1500 0.9750 /
0.1800 514.7 1.2070 0.9100 /
0.3710 1014.7 1.2950 0.8300 /
0.6360 2014.7 1.4350 0.6950 /
0.7750 2514.7 1.5000 0.6410 /
0.9300 3014.7 1.5650 0.5940 /
1.2700 4014.7 1.6950 0.5100
9014.7 1.5790 0.7400 /
1.6180 5014.7 1.8270 0.4490
9014.7 1.7370 0.6310 /
/
SOLUTION
SWAT
100*0.0 /
SGAS
100*0.0 /
PRESSURE
100*300.0 /
SUMMARY
SCHEDULE

View File

@ -0,0 +1,87 @@
-- Most of the following sections are not actually needed by the test,
-- but it is required by the Eclipse reference manual that they are
-- present. Also, the higher level opm-parser classes
-- (i.e. Opm::EclipseState et al.) assume that they are present.
-------------------------------------
RUNSPEC
WATER
OIL
GAS
DIMENS
1 1 20 /
TABDIMS
1 1 40 20 1 20 /
EQLDIMS
-- NTEQUL
1 /
-------------------------------------
GRID
-- Opm::EclipseState assumes that _some_ grid gets defined, so let's
-- specify a fake one...
DXV
1*1 /
DYV
1*1 /
DZV
20*5 /
DEPTHZ
4*0.0 /
PORO
20*0.3 /
PERMX
20*500 /
-------------------------------------
PROPS
PVDO
100 1.0 1.0
200 0.9 1.0
/
PVDG
100 0.010 0.1
200 0.005 0.2
/
PVTW
1.0 1.0 4.0E-5 0.96 0.0
/
SWOF
0.2 0 1 0.4
1 1 0 0.1
/
SGOF
0 0 1 0.2
0.8 1 0 0.5
/
DENSITY
700 1000 1
/
-------------------------------------
SOLUTION
EQUIL
50 150 50 0.25 20 0.35 1* 1* 0
/
-------------------------------------
SCHEDULE
-- empty section

View File

@ -0,0 +1,128 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 20
/
WELLDIMS
30 10 2 30 /
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DXV
1.0
/
DYV
1.0
/
DZV
20*5.0
/
PORO
20*0.2
/
PERMZ
20*1.0
/
PERMY
20*100.0
/
PERMX
20*100.0
/
BOX
1 1 1 1 1 1 /
TOPS
0.0
/
PROPS ======
PVDO
100 1.0 1.0
200 0.9 1.0
/
PVDG
100 0.010 0.1
200 0.005 0.2
/
SWOF
0.2 0 1 0.9
1 1 0 0.1
/
SGOF
0 0 1 0.2
0.8 1 0 0.5
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1* 1* 0
/
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RESTART=3' 'NEWTON=2' /
END

View File

@ -0,0 +1,102 @@
-- Most of the following sections are not actually needed by the test,
-- but it is required by the Eclipse reference manual that they are
-- present. Also, the higher level opm-parser classes
-- (i.e. Opm::EclipseState et al.) assume that they are present.
-------------------------------------
RUNSPEC
WATER
OIL
GAS
DIMENS
1 1 20 /
TABDIMS
1 1 40 20 1 20 /
EQLDIMS
-- NTEQUL
1 /
START
1 'JAN' 2015 /
-------------------------------------
GRID
-- Opm::EclipseState assumes that _some_ grid gets defined, so let's
-- specify a fake one...
DXV
1 /
DYV
1 /
DZ
20*5 /
TOPS
0 /
PORO
20*0.3 /
PERMX
20*500 /
PERMZ
20*500 /
-------------------------------------
PROPS
ROCK
14.7 3E-6 /
PVDO
100 1.0 1.0
200 0.9 1.0
/
PVDG
100 0.010 0.1
200 0.005 0.2
/
PVTW
1.0 1.0 4.0E-5 0.96 0.0
/
SWOF
0.2 0 1 0.4
1 1 0 0.1
/
SGOF
0 0 1 0.2
0.8 1 0 0.5
/
DENSITY
700 1000 1
/
SWATINIT
5*0
10*0.5
5*1 /
-------------------------------------
SOLUTION
EQUIL
50 150 50 0.25 20 0.35 1* 1* 0
/
RPTSOL
SWATINIT SWAT SGAS SOIL PCOG PCOW
/
-------------------------------------
SCHEDULE
-- empty section

View File

@ -0,0 +1,88 @@
-- Most of the following sections are not actually needed by the test,
-- but it is required by the Eclipse reference manual that they are
-- present. Also, the higher level opm-parser classes
-- (i.e. Opm::EclipseState et al.) assume that they are present.
-------------------------------------
RUNSPEC
WATER
GAS
OIL
METRIC
DIMENS
1 1 10 /
TABDIMS
1 1 40 20 1 20 /
EQLDIMS
-- NTEQUL
1 /
GRID
-- Opm::EclipseState assumes that _some_ grid gets defined, so let's
-- specify a fake one...
DXV
1*1 /
DYV
1*1 /
DZV
10*1 /
TOPS
1*0.0 /
PORO
10*0.3 /
PERMX
10*500 /
-------------------------------------
PROPS
PVDO
100 1.0 1.0
200 0.5 1.0
/
PVDG
100 0.05 0.1
200 0.02 0.2
/
PVTW
1.0 1.0 4.0E-5 0.96 0.0
/
SWOF
0 0 1 0
1 1 0 0
/
SGOF
0 0 1 0
1 1 0 0
/
DENSITY
700 1000 10
/
-------------------------------------
SOLUTION
EQUIL
5 150 5 0 2 0 1* 1* 0
/
-------------------------------------
SCHEDULE
-- empty section

View File

@ -0,0 +1,152 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
DISGAS
VAPOIL
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 20
/
WELLDIMS
30 10 2 30 /
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DX
20*1.0
/
DY
20*1.0
/
DZ
20*5.0
/
TOPS
0.0
/
PORO
20*0.2 /
PERMZ
20*1.0
/
PERMY
20*100.0
/
PERMX
20*100.0
/
BOX
1 1 1 1 1 1 /
PROPS ======
PVTO
-- Rs Pbub Bo Vo
0 1. 1.0000 1.20 /
20 40. 1.0120 1.17 /
40 80. 1.0255 1.14 /
60 120. 1.0380 1.11 /
80 160. 1.0510 1.08 /
100 200. 1.0630 1.06 /
120 240. 1.0750 1.03 /
140 280. 1.0870 1.00 /
160 320. 1.0985 .98 /
180 360. 1.1100 .95 /
200 400. 1.1200 .94
500. 1.1189 .94 /
/
PVTG
-- Pg Rv Bg Vg
100 0.0001 0.010 0.1
0.0 0.0104 0.1 /
200 0.0004 0.005 0.2
0.0 0.0054 0.2 /
/
SWOF
0.2 0 1 0.9
1 1 0 0.1
/
SGOF
0 0 1 0.2
0.8 1 0 0.5
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1 1 0
/
PBVD
0 1.0
50 150. /
PDVD
50. 150.
100. 100 /
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
END

258
tests/test_ecl_output.cc Normal file
View File

@ -0,0 +1,258 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#include "config.h"
#include <ebos/equil/equilibrationhelpers.hh>
#include <ebos/eclproblem.hh>
#include <ewoms/common/start.hh>
#include <opm/grid/UnstructuredGrid.h>
#include <opm/grid/GridManager.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/output/eclipse/Summary.hpp>
#include <ebos/collecttoiorank.hh>
#include <ebos/ecloutputblackoilmodule.hh>
#include <ebos/eclwriter.hh>
#if HAVE_DUNE_FEM
#include <dune/fem/misc/mpimanager.hh>
#else
#include <dune/common/parallel/mpihelper.hh>
#endif
#include <array>
#include <iostream>
#include <limits>
#include <memory>
#include <numeric>
#include <sstream>
#include <string>
#include <vector>
#include <string.h>
#define CHECK(value, expected) \
{ \
if ((value) != (expected)) \
std::abort(); \
}
#define CHECK_CLOSE(value, expected, reltol) \
{ \
if (std::fabs((expected) - (value)) > 1e-14 && \
std::fabs(((expected) - (value))/((expected) + (value))) > reltol) \
{ \
std::cout << "Test failure: "; \
std::cout << "expected value " << expected << " is not close to value " << value << std::endl; \
std::abort(); \
} \
} \
#define REQUIRE(cond) \
{ \
if (!(cond)) \
std::abort(); \
}
BEGIN_PROPERTIES
NEW_TYPE_TAG(TestEclOutputTypeTag, INHERITS_FROM(BlackOilModel, EclBaseProblem));
SET_BOOL_PROP(TestEclOutputTypeTag, EnableGravity, false);
SET_BOOL_PROP(TestEclOutputTypeTag, EnableAsyncEclOutput, false);
END_PROPERTIES
static const int day = 24 * 60 * 60;
template <class TypeTag>
std::unique_ptr<typename GET_PROP_TYPE(TypeTag, Simulator)>
initSimulator(const char *filename)
{
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
std::string filenameArg = "--ecl-deck-file-name=";
filenameArg += filename;
const char* argv[] = {
"test_equil",
filenameArg.c_str()
};
Ewoms::setupParameters_<TypeTag>(/*argc=*/sizeof(argv)/sizeof(argv[0]), argv, /*registerParams=*/false);
return std::unique_ptr<Simulator>(new Simulator);
}
ERT::ert_unique_ptr<ecl_sum_type, ecl_sum_free> readsum(const std::string& base);
ERT::ert_unique_ptr<ecl_sum_type, ecl_sum_free> readsum(const std::string& base)
{
return ERT::ert_unique_ptr<ecl_sum_type, ecl_sum_free>(
ecl_sum_fread_alloc_case(base.c_str(), ":"));
}
void test_summary();
void test_summary()
{
typedef typename TTAG(TestEclOutputTypeTag) TypeTag;
const std::string filename = "SUMMARY_DECK_NON_CONSTANT_POROSITY.DATA";
const std::string casename = "SUMMARY_DECK_NON_CONSTANT_POROSITY";
auto simulator = initSimulator<TypeTag>(filename.data());
typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard;
typedef Ewoms::CollectDataToIORank< Vanguard > CollectDataToIORankType;
CollectDataToIORankType collectToIORank(simulator->vanguard());
Ewoms::EclOutputBlackOilModule<TypeTag> eclOutputModule(*simulator, collectToIORank);
typedef Ewoms::EclWriter<TypeTag> EclWriterType;
// create the actual ECL writer
std::unique_ptr<EclWriterType> eclWriter = std::unique_ptr<EclWriterType>(new EclWriterType(*simulator));
simulator->model().applyInitialSolution();
Opm::data::Wells dw;
bool substep = false;
simulator->setEpisodeIndex(0);
eclWriter->writeOutput(substep);
simulator->setEpisodeIndex(1);
eclWriter->writeOutput(substep);
simulator->setEpisodeIndex(2);
eclWriter->writeOutput(substep);
auto res = readsum( casename );
const auto* resp = res.get();
// fpr = sum_ (p * hcpv ) / hcpv, hcpv = pv * (1 - sw)
const double fpr = ( (3 * 0.1 + 8 * 0.2) * 500 * (1 - 0.2) ) / ( (500*0.1 + 500*0.2) * (1 - 0.2));
CHECK_CLOSE( fpr, ecl_sum_get_field_var( resp, 1, "FPR" ) , 1e-5 );
// foip = sum_ (b * s * pv), rs == 0;
const double foip = ( (0.3 * 0.1 + 0.8 * 0.2) * 500 * (1 - 0.2) );
CHECK_CLOSE(foip, ecl_sum_get_field_var( resp, 1, "FOIP" ), 1e-3 );
// fgip = sum_ (b * pv * s), sg == 0;
const double fgip = 0.0;
CHECK_CLOSE(fgip, ecl_sum_get_field_var( resp, 1, "FGIP" ), 1e-3 );
// fgip = sum_ (b * pv * s),
const double fwip = 1.0/1000 * ( 0.1 + 0.2) * 500 * 0.2;
CHECK_CLOSE(fwip, ecl_sum_get_field_var( resp, 1, "FWIP" ), 1e-3 );
// region 1
// rpr = sum_ (p * hcpv ) / hcpv, hcpv = pv * (1 - sw)
const double rpr1 = ( 2.5 * 0.1 * 400 * (1 - 0.2) ) / (400*0.1 * (1 - 0.2));
CHECK_CLOSE( rpr1, ecl_sum_get_general_var( resp, 1, "RPR:1" ) , 1e-5 );
// roip = sum_ (b * s * pv) // rs == 0;
const double roip1 = ( 0.25 * 0.1 * 400 * (1 - 0.2) );
CHECK_CLOSE(roip1, ecl_sum_get_general_var( resp, 1, "ROIP:1" ), 1e-3 );
// region 2
// rpr = sum_ (p * hcpv ) / hcpv, hcpv = pv * (1 - sw)
const double rpr2 = ( (5 * 0.1 * 100 + 6 * 0.2 * 100) * (1 - 0.2) ) / ( (100*0.1 + 100*0.2) * (1 - 0.2));
CHECK_CLOSE( rpr2, ecl_sum_get_general_var( resp, 1, "RPR:2" ) , 1e-5 );
// roip = sum_ (b * s * pv) // rs == 0;
const double roip2 = ( (0.5 * 0.1 * 100 + 0.6 * 0.2 * 100) * (1 - 0.2) );
CHECK_CLOSE(roip2, ecl_sum_get_general_var( resp, 1, "ROIP:2" ), 1e-3 );
}
void test_readWriteWells();
void test_readWriteWells()
{
using opt = Opm::data::Rates::opt;
Opm::data::Rates r1, r2, rc1, rc2, rc3;
r1.set( opt::wat, 5.67 );
r1.set( opt::oil, 6.78 );
r1.set( opt::gas, 7.89 );
r2.set( opt::wat, 8.90 );
r2.set( opt::oil, 9.01 );
r2.set( opt::gas, 10.12 );
rc1.set( opt::wat, 20.41 );
rc1.set( opt::oil, 21.19 );
rc1.set( opt::gas, 22.41 );
rc2.set( opt::wat, 23.19 );
rc2.set( opt::oil, 24.41 );
rc2.set( opt::gas, 25.19 );
rc3.set( opt::wat, 26.41 );
rc3.set( opt::oil, 27.19 );
rc3.set( opt::gas, 28.41 );
Opm::data::Well w1, w2;
w1.rates = r1;
w1.bhp = 1.23;
w1.temperature = 3.45;
w1.control = 1;
/*
* the connection keys (active indices) and well names correspond to the
* input deck. All other entries in the well structures are arbitrary.
*/
w1.connections.push_back( { 88, rc1, 30.45, 123.45, 0.0, 0.0, 0.0 } );
w1.connections.push_back( { 288, rc2, 33.19, 67.89, 0.0, 0.0, 0.0 } );
w2.rates = r2;
w2.bhp = 2.34;
w2.temperature = 4.56;
w2.control = 2;
w2.connections.push_back( { 188, rc3, 36.22, 19.28, 0.0, 0.0, 0.0 } );
Opm::data::Wells wellRates;
wellRates["OP_1"] = w1;
wellRates["OP_2"] = w2;
typedef Dune :: Point2PointCommunicator< Dune :: SimpleMessageBuffer > P2PCommunicatorType;
typedef typename P2PCommunicatorType :: MessageBufferType MessageBufferType;
MessageBufferType buffer;
wellRates.write(buffer);
Opm::data::Wells wellRatesCopy;
wellRatesCopy.read(buffer);
CHECK( wellRatesCopy.get( "OP_1" , opt::wat) , wellRates.get( "OP_1" , opt::wat));
CHECK( wellRatesCopy.get( "OP_2" , 188 , opt::wat) , wellRates.get( "OP_2" , 188 , opt::wat));
}
int main(int argc, char** argv)
{
#if HAVE_DUNE_FEM
Dune::Fem::MPIManager::initialize(argc, argv);
#else
Dune::MPIHelper::instance(argc, argv);
#endif
typedef TTAG(TestEclOutputTypeTag) TypeTag;
Ewoms::registerAllParameters_<TypeTag>();
test_summary();
test_readWriteWells();
return 0;
}

1127
tests/test_equil.cc Normal file

File diff suppressed because it is too large Load Diff