EclSimulator: make compile with dune-alugrid. For that reason the Cartesian cell id has

been moved to the grid manager class.
This commit is contained in:
Robert K 2015-01-19 16:22:34 +01:00
parent 83f335077a
commit 779e7736d1
3 changed files with 47 additions and 38 deletions

View File

@ -503,12 +503,12 @@ public:
if (!deck->hasKeyword("PVTNUM"))
return 0;
const auto &grid = this->simulator().gridManager().grid();
const auto &cartesianCellId = this->simulator().gridManager().cartesianCellId();
// this is quite specific to the ECFV discretization. But so is everything in an
// ECL deck, i.e., we don't need to care here...
int compressedDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
int cartesianDofIdx = grid.globalCell()[compressedDofIdx];
int cartesianDofIdx = cartesianCellId[compressedDofIdx];
return deck->getKeyword("PVTNUM")->getIntData()[cartesianDofIdx] - 1;
}
@ -638,7 +638,7 @@ private:
{
auto deck = this->simulator().gridManager().deck();
auto eclState = this->simulator().gridManager().eclState();
const auto &grid = this->simulator().gridManager().grid();
const auto &cartesianCellId = this->simulator().gridManager().cartesianCellId();
size_t numDof = this->model().numGridDof();
@ -663,7 +663,7 @@ private:
permzData = eclState->getDoubleGridProperty("PERMZ")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = grid.globalCell()[dofIdx];
int cartesianElemIdx = cartesianCellId[dofIdx];
intrinsicPermeability_[dofIdx] = 0.0;
intrinsicPermeability_[dofIdx][0][0] = permxData[cartesianElemIdx];
intrinsicPermeability_[dofIdx][1][1] = permyData[cartesianElemIdx];
@ -686,7 +686,7 @@ private:
eclState->getDoubleGridProperty("PORO")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = grid.globalCell()[dofIdx];
int cartesianElemIdx = cartesianCellId[dofIdx];
porosity_[dofIdx] = poroData[cartesianElemIdx];
}
}
@ -701,7 +701,7 @@ private:
eclState->getDoubleGridProperty("NTG")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = grid.globalCell()[dofIdx];
int cartesianElemIdx = cartesianCellId[dofIdx];
porosity_[dofIdx] *= ntgData[cartesianElemIdx];
}
}
@ -712,7 +712,7 @@ private:
eclState->getDoubleGridProperty("MULTPV")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = grid.globalCell()[dofIdx];
int cartesianElemIdx = cartesianCellId[dofIdx];
porosity_[dofIdx] *= multpvData[cartesianElemIdx];
}
}
@ -776,7 +776,7 @@ private:
materialParamTableIdx_.resize(numDof);
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = grid.globalCell()[dofIdx];
int cartesianElemIdx = cartesianCellId[dofIdx];
// make sure that all values are in the correct range
assert(1 <= satnumData[dofIdx]);
@ -822,7 +822,7 @@ private:
void readInitialCondition_()
{
const auto deck = this->simulator().gridManager().deck();
const auto &grid = this->simulator().gridManager().grid();
const auto &cartesianCellId = this->simulator().gridManager().cartesianCellId();
if (!deck->hasKeyword("SWAT") ||
!deck->hasKeyword("SGAS"))
@ -864,7 +864,7 @@ private:
// make sure that the size of the data arrays is correct
#ifndef NDEBUG
const auto &cartSize = grid.logicalCartesianSize();
const auto &cartSize = this->simulator().gridManager().logicalCartesianSize();
size_t numCartesianCells = cartSize[0] * cartSize[1] * cartSize[2];
assert(waterSaturationData.size() == numCartesianCells);
assert(gasSaturationData.size() == numCartesianCells);
@ -876,7 +876,7 @@ private:
for (size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
auto &dofFluidState = initialFluidStates_[dofIdx];
size_t cartesianDofIdx = grid.globalCell()[dofIdx];
size_t cartesianDofIdx = cartesianCellId[dofIdx];
assert(0 <= cartesianDofIdx);
assert(cartesianDofIdx <= numCartesianCells);
@ -928,7 +928,7 @@ private:
if (RsReal > RsSat) {
std::array<int, 3> ijk;
grid.getIJK(dofIdx, ijk);
// grid.getIJK(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"

View File

@ -97,14 +97,15 @@ public:
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& entity = *elemIt;
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
int elemIdx = elementMapper.index(elemIt);
int elemIdx = elementMapper.index( entity );
#else
int elemIdx = elementMapper.map(elemIt);
int elemIdx = elementMapper.map( entity );
#endif
// get the geometry of the current element
const auto& geom = elemIt->geometry();
const auto& geom = entity.geometry();
// compute the axis specific "centroids" used for the
// transmissibilities
@ -153,19 +154,25 @@ public:
// compute the transmissibilities for all intersections
elemIt = gridView.template begin</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
auto isIt = elemIt->ileafbegin();
const auto& isEndIt = elemIt->ileafend();
const auto& entity = *elemIt;
auto isIt = gridView.ibegin( entity );
const auto& isEndIt = gridView.iend( entity );
for (; isIt != isEndIt; ++ isIt) {
// store intersection, this might be costly
const auto& intersection = *isIt;
// ignore boundary intersections for now (TODO?)
if (isIt->boundary())
if (intersection.boundary())
continue;
const auto& inside = intersection.inside();
const auto& outside = intersection.outside();
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
int insideElemIdx = elementMapper.index(*isIt->inside());
int outsideElemIdx = elementMapper.index(*isIt->outside());
int insideElemIdx = elementMapper.index( *inside );
int outsideElemIdx = elementMapper.index( *outside);
#else
int insideElemIdx = elementMapper.map(*isIt->inside());
int outsideElemIdx = elementMapper.map(*isIt->outside());
int insideElemIdx = elementMapper.map( *inside );
int outsideElemIdx = elementMapper.map( *outside);
#endif
// we only need to calculate a face's transmissibility
@ -175,25 +182,25 @@ public:
// local indices of the faces of the inside and
// outside elements which contain the intersection
int insideFaceIdx = isIt->indexInInside();
int outsideFaceIdx = isIt->indexInOutside();
int insideFaceIdx = intersection.indexInInside();
int outsideFaceIdx = intersection.indexInOutside();
Scalar halfTrans1;
Scalar halfTrans2;
computeHalfTrans_(halfTrans1,
*isIt,
intersection,
insideFaceIdx,
distanceVector_(*isIt,
isIt->indexInInside(),
distanceVector_(intersection,
intersection.indexInInside(),
insideElemIdx,
axisCentroids),
problem.intrinsicPermeability(insideElemIdx));
computeHalfTrans_(halfTrans2,
*isIt,
intersection,
outsideFaceIdx,
distanceVector_(*isIt,
isIt->indexInOutside(),
distanceVector_(intersection,
intersection.indexInOutside(),
outsideElemIdx,
axisCentroids),
problem.intrinsicPermeability(outsideElemIdx));

View File

@ -539,6 +539,7 @@ protected:
void updateWellTopology_(int reportStepIdx, const WellCompletionsMap& wellCompletions)
{
auto& model = simulator_.model();
const auto& cartesianCellId = simulator_.gridManager().cartesianCellId();
// first, remove all wells from the reservoir
model.clearAuxiliaryModules();
@ -548,20 +549,20 @@ protected:
(*wellIt)->clear();
// tell the active wells which DOFs they contain
const auto& grid = simulator_.gridManager().grid();
const auto gridView = simulator_.gridManager().gridView();
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) {
if (elemIt->partitionType() != Dune::InteriorEntity)
const auto& entity = *elemIt;
if (entity.partitionType() != Dune::InteriorEntity)
continue; // non-local entities need to be skipped
elemCtx.updateStencil(elemIt);
elemCtx.updateStencil( entity );
for (int dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx) {
int globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
int cartesianDofIdx = grid.globalCell()[globalDofIdx];
int cartesianDofIdx = cartesianCellId[globalDofIdx];
if (wellCompletions.count(cartesianDofIdx) == 0)
// the current DOF is not contained in any well, so we must skip
@ -626,21 +627,22 @@ protected:
{
// associate the well completions with grid cells and register them in the
// Peaceman well object
const Grid& grid = simulator_.gridManager().grid();
const GridView gridView = simulator_.gridManager().gridView();
const auto& cartesianCellId = simulator_.gridManager().cartesianCellId();
ElementContext elemCtx(simulator_);
auto elemIt = gridView.template begin</*codim=*/0>();
const auto elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
if (elemIt->partitionType() != Dune::InteriorEntity)
const auto& entity = *elemIt;
if (entity.partitionType() != Dune::InteriorEntity)
continue; // non-local entities need to be skipped
elemCtx.updateStencil(elemIt);
elemCtx.updateStencil( entity );
for (int dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx) {
int globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
int cartesianDofIdx = grid.globalCell()[globalDofIdx];
int cartesianDofIdx = cartesianCellId[globalDofIdx];
if (wellCompletions.count(cartesianDofIdx) == 0)
// the current DOF is not contained in any well, so we must skip