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

CartesianIndexMapper: use the CartesianIndexMapper to applied the mapping from the

Browse Source 
Cartesian coordinate to the flat index and vice versa.
This commit is contained in:
Robert Kloefkorn
2015-09-10 15:55:54 +02:00
parent 58aeddf431
commit f6b3bbc386
5 changed files with 249 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/ebos/eclequilinitializer.hh
View File

@@ -135,7 +135,7 @@ public:
// set the phase pressures. the Opm::BlackoilState only provides the oil
// phase pressure, so we need to calculate the other phases' pressures
// ourselfs.
Scalar pC[numPhases];
Dune::FieldVector< Scalar, numPhases > pC( 0 );
const auto& matParams = simulator.problem().materialLawParams(elemIdx);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
Scalar po = opmBlackoilState.pressure()[elemIdx];

27
applications/ebos/eclproblem.hh
View File

@@ -596,8 +596,7 @@ public:
const auto& gridManager = this->simulator().gridManager();
int cartesianDofIdx = gridManager.cartesianCellId(elemIdx);
int cartesianDofIdx = gridManager.cartesianIndex(elemIdx);
return deck->getKeyword("PVTNUM")->getIntData()[cartesianDofIdx] - 1;
}
@@ -716,7 +715,7 @@ private:
eclState->getIntGridProperty("PVTNUM")->getData();
rockTableIdx_.resize(gridManager.gridView().size(/*codim=*/0));
for (size_t elemIdx = 0; elemIdx < rockTableIdx_.size(); ++ elemIdx) {
int cartElemIdx = gridManager.cartesianCellId(elemIdx);
int cartElemIdx = gridManager.cartesianIndex(elemIdx);
// reminder: Eclipse uses FORTRAN-style indices
rockTableIdx_[elemIdx] = pvtnumData[cartElemIdx] - 1;
@@ -751,7 +750,7 @@ private:
permzData = eclState->getDoubleGridProperty("PERMZ")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = gridManager.cartesianCellId(dofIdx);
int cartesianElemIdx = gridManager.cartesianIndex(dofIdx);
intrinsicPermeability_[dofIdx] = 0.0;
intrinsicPermeability_[dofIdx][0][0] = permxData[cartesianElemIdx];
intrinsicPermeability_[dofIdx][1][1] = permyData[cartesianElemIdx];
@@ -779,7 +778,7 @@ private:
eclState->getDoubleGridProperty("PORO")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = gridManager.cartesianCellId(dofIdx);
int cartesianElemIdx = gridManager.cartesianIndex(dofIdx);
porosity_[dofIdx] = poroData[cartesianElemIdx];
}
}
@@ -791,7 +790,7 @@ private:
eclState->getDoubleGridProperty("PORV")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = gridManager.cartesianCellId(dofIdx);
int cartesianElemIdx = gridManager.cartesianIndex(dofIdx);
if (std::isfinite(porvData[cartesianElemIdx])) {
Scalar dofVolume = this->simulator().model().dofTotalVolume(dofIdx);
porosity_[dofIdx] = porvData[cartesianElemIdx]/dofVolume;
@@ -805,7 +804,7 @@ private:
eclState->getDoubleGridProperty("NTG")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx)
porosity_[dofIdx] *= ntgData[gridManager.cartesianCellId(dofIdx)];
porosity_[dofIdx] *= ntgData[gridManager.cartesianIndex(dofIdx)];
}
// apply the MULTPV keyword to the porosity
@@ -814,13 +813,13 @@ private:
eclState->getDoubleGridProperty("MULTPV")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx)
porosity_[dofIdx] *= multpvData[gridManager.cartesianCellId(dofIdx)];
porosity_[dofIdx] *= multpvData[gridManager.cartesianIndex(dofIdx)];
}
// the fluid-matrix interactions for ECL problems are dealt with by a separate class
std::vector<int> compressedToCartesianElemIdx(numDof);
for (unsigned elemIdx = 0; elemIdx < numDof; ++elemIdx)
compressedToCartesianElemIdx[elemIdx] = gridManager.cartesianCellId(elemIdx);
compressedToCartesianElemIdx[elemIdx] = gridManager.cartesianIndex(elemIdx);
materialLawManager_ = std::make_shared<EclMaterialLawManager>();
materialLawManager_->initFromDeck(deck, eclState, compressedToCartesianElemIdx);
@@ -1042,7 +1041,7 @@ private:
// make sure that the size of the data arrays is correct
#ifndef NDEBUG
const auto &cartSize = this->simulator().gridManager().logicalCartesianSize();
const auto &cartSize = this->simulator().gridManager().cartesianDimensions();
size_t numCartesianCells = cartSize[0] * cartSize[1] * cartSize[2];
assert(waterSaturationData.size() == numCartesianCells);
assert(gasSaturationData.size() == numCartesianCells);
@@ -1057,7 +1056,7 @@ private:
for (size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
auto &dofFluidState = initialFluidStates_[dofIdx];
size_t cartesianDofIdx = gridManager.cartesianCellId(dofIdx);
size_t cartesianDofIdx = gridManager.cartesianIndex(dofIdx);
assert(0 <= cartesianDofIdx);
assert(cartesianDofIdx <= numCartesianCells);
@@ -1088,7 +1087,7 @@ private:
// this assumes that capillary pressures only depend on the phase saturations
// and possibly on temperature. (this is always the case for ECL problems.)
Scalar pc[numPhases];
Dune::FieldVector< Scalar, numPhases > pc( 0 );
const auto& matParams = materialLawParams(dofIdx);
MaterialLaw::capillaryPressures(pc, matParams, dofFluidState);
Valgrind::CheckDefined(oilPressure);
@@ -1127,7 +1126,7 @@ private:
if (RsReal > RsSat) {
std::array<int, 3> ijk;
gridManager.getIJK(dofIdx, ijk);
gridManager.cartesianCoordinate(dofIdx, ijk);
std::cerr << "Warning: The specified amount gas (R_s = " << RsReal << ") is more"
<< " than the maximium\n"
<< " amount which can be dissolved in oil"
@@ -1167,7 +1166,7 @@ private:
if (RvReal > RvSat) {
std::array<int, 3> ijk;
gridManager.getIJK(dofIdx, ijk);
gridManager.cartesianCoordinate(dofIdx, ijk);
std::cerr << "Warning: The specified amount oil (R_v = " << RvReal << ") is more"
<< " than the maximium\n"
<< " amount which can be dissolved in gas"

4
applications/ebos/ecltransmissibility.hh
View File

@@ -184,8 +184,8 @@ public:
if (insideElemIdx > outsideElemIdx)
continue;
int cartesianElemIdxInside = gridManager.cartesianCellId(insideElemIdx);
int cartesianElemIdxOutside = gridManager.cartesianCellId(outsideElemIdx);
int cartesianElemIdxInside = gridManager.cartesianIndex(insideElemIdx);
int cartesianElemIdxOutside = gridManager.cartesianIndex(outsideElemIdx);
// local indices of the faces of the inside and
// outside elements which contain the intersection

17
applications/ebos/eclwellmanager.hh
View File

@@ -656,7 +656,7 @@ protected:
elemCtx.updateStencil(elem);
for (int dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx) {
int globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
int cartesianDofIdx = gridManager.cartesianCellId(globalDofIdx);
int cartesianDofIdx = gridManager.cartesianIndex(globalDofIdx);
if (wellCompletions.count(cartesianDofIdx) == 0)
// the current DOF is not contained in any well, so we must skip
@@ -683,9 +683,8 @@ protected:
auto deckSchedule = eclStatePtr->getSchedule();
auto eclGrid = eclStatePtr->getEclipseGrid();
int nx = eclGrid->getNX();
int ny = eclGrid->getNY();
//int nz = eclGrid->getNZ();
assert( eclGrid->getNX() == simulator_.gridManager().cartesianDimensions()[ 0 ] );
assert( eclGrid->getNY() == simulator_.gridManager().cartesianDimensions()[ 1 ] );
// compute the mapping from logically Cartesian indices to the well the
// respective completion.
@@ -701,11 +700,17 @@ protected:
continue;
}
std::array<int, 3> cartesianCoordinate;
// set the well parameters defined by the current set of completions
Opm::CompletionSetConstPtr completionSet = deckWell->getCompletions(reportStepIdx);
for (size_t complIdx = 0; complIdx < completionSet->size(); complIdx ++) {
Opm::CompletionConstPtr completion = completionSet->get(complIdx);
int cartIdx = completion->getI() + completion->getJ()*nx + completion->getK()*nx*ny;
cartesianCoordinate[ 0 ] = completion->getI();
cartesianCoordinate[ 1 ] = completion->getJ();
cartesianCoordinate[ 2 ] = completion->getK();
const int cartIdx = simulator_.gridManager().cartesianIndex( cartesianCoordinate );
assert( cartIdx == (completion->getI() + completion->getJ()*eclGrid->getNX()
+ completion->getK()*eclGrid->getNX()*eclGrid->getNY() ) );
// in this code we only support each cell to be part of at most a single
// well. TODO (?) change this?
@@ -748,7 +753,7 @@ protected:
elemCtx.updateStencil(elem);
for (int dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx) {
int globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
int cartesianDofIdx = gridManager.cartesianCellId(globalDofIdx);
int cartesianDofIdx = gridManager.cartesianIndex(globalDofIdx);
if (wellCompletions.count(cartesianDofIdx) == 0)
// the current DOF is not contained in any well, so we must skip

222
ewoms/common/cartesianindexmapper.hh Normal file
View File

@@ -0,0 +1,222 @@
#ifndef OPM_CARTESIANINDEXMAPPER_HEADER
#define OPM_CARTESIANINDEXMAPPER_HEADER
#include <array>
#include <cassert>
#include <dune/common/exceptions.hh>
#include <dune/grid/common/datahandleif.hh>
#include <dune/grid/utility/persistentcontainer.hh>
namespace Dune
{
/** \brief Interface class to access the logical Cartesian grid as used in industry
standard simulator decks.
*/
template< class Grid >
class CartesianIndexMapper
{
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
{
int idx_;
public:
GlobalCellIndex() : idx_(-1) {}
GlobalCellIndex& operator= ( const int index ) { idx_ = index; return *this; }
int index() const { return idx_; }
};
typedef typename Dune::PersistentContainer< Grid, GlobalCellIndex > GlobalIndexContainer;
GridView gridView_;
GlobalIndexContainer globalIndex_;
std::vector<int>& cartesianIndex_;
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();
//std::cout << "CartesianIndex " << cartesianIndex_.size() << std::endl;
//for( size_t i=0; i<cartesianIndex_.size(); ++i )
// std::cout << "cart[ " << i << " ] = " << cartesianIndex_[ i ] << std::endl;
}
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();
for(auto it = gridView_.template begin<0>(),
end = gridView_.template end<0>(); it != end; ++it, ++idx )
{
globalIndex_[ *it ] = *idx;
}
}
// update vector from given persistent container
void finalize()
{
std::vector< int > newIndex ;
newIndex.reserve( gridView_.indexSet().size( 0 ) );
for(auto it = gridView_.template begin<0>(),
end = gridView_.template end<0>(); it != end; ++it)
{
newIndex.push_back( globalIndex_[ *it ].index() ) ;
}
cartesianIndex_.swap( newIndex );
}
};
public:
/** \brief dimension of the grid */
static const int dimension = Grid :: dimension ;
/** \brief constructor taking grid */
CartesianIndexMapper( 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 >=2 )
{
for( int d=0; d<dimension-2; ++d )
{
coords[d] = gc % cartesianDimensions()[d]; gc /= cartesianDimensions()[d];
}
coords[dimension-2] = gc % cartesianDimensions()[dimension-2];
coords[dimension-1] = gc / cartesianDimensions()[dimension-1];
}
else
coords[ 0 ] = gc ;
}
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 ;
}
protected:
const Grid& grid_;
const std::array<int, dimension> cartesianDimensions_;
std::vector<int> cartesianIndex_;
const int cartesianSize_ ;
};
} // end namespace Opm
#endif