From 14418092e13dd7288fcaf598c0bd55d40369b9ee Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 6 Sep 2016 15:16:03 +0200 Subject: [PATCH] ebos: make the porosity handling more general i.e., it now uses the PORV grid property from opm-parser and does the accumulation of the disabled cells manually. This patch should be equivalent to the opm-simulators PR #806 ( https://github.com/OPM/opm-simulators/pull/806 ). --- applications/ebos/eclproblem.hh | 91 ++++++++++++++++----------------- 1 file changed, 45 insertions(+), 46 deletions(-) diff --git a/applications/ebos/eclproblem.hh b/applications/ebos/eclproblem.hh index 9149c3347..a90fbc19e 100644 --- a/applications/ebos/eclproblem.hh +++ b/applications/ebos/eclproblem.hh @@ -812,15 +812,15 @@ private: const auto& gridManager = this->simulator().gridManager(); Opm::DeckConstPtr deck = gridManager.deck(); Opm::EclipseStateConstPtr eclState = gridManager.eclState(); + Opm::EclipseGridConstPtr eclGrid = eclState->getInputGrid(); const auto& props = eclState->get3DProperties(); size_t numDof = this->model().numGridDof(); - intrinsicPermeability_.resize(numDof); - porosity_.resize(numDof); - //////////////////////////////// - // permeability + // intrinsic permeability + + intrinsicPermeability_.resize(numDof); // read the intrinsic permeabilities from the eclState. Note that all arrays // provided by eclState are one-per-cell of "uncompressed" grid, whereas the @@ -853,62 +853,61 @@ private: //////////////////////////////// - // compute the porosity - if (!props.hasDeckDoubleGridProperty("PORO") && !props.hasDeckDoubleGridProperty("PORV")) - OPM_THROW(std::runtime_error, - "Can't read the porosity from the ECL state object. " - "(The PORO and PORV keywords are missing)"); + // porosity + porosity_.resize(numDof); - if (props.hasDeckDoubleGridProperty("PORO")) { - const std::vector &poroData = - props.getDoubleGridProperty("PORO").getData(); + const std::vector &porvData = + props.getDoubleGridProperty("PORV").getData(); + const std::vector &actnumData = + props.getIntGridProperty("ACTNUM").getData(); - for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) { - unsigned cartesianElemIdx = gridManager.cartesianIndex(dofIdx); - porosity_[dofIdx] = poroData[cartesianElemIdx]; - } - } + int nx = eclGrid->getNX(); + int ny = eclGrid->getNY(); + for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) { + unsigned cartElemIdx = gridManager.cartesianIndex(dofIdx); + Scalar poreVolume = porvData[cartElemIdx]; - // apply the NTG keyword to the porosity - if (props.hasDeckDoubleGridProperty("NTG")) { - const std::vector &ntgData = - props.getDoubleGridProperty("NTG").getData(); + // sum up the pore volume of the active cell and all inactive ones above it + // which were disabled due to their pore volume being too small + if (eclGrid->getMinpvMode() == Opm::MinpvMode::ModeEnum::OpmFIL) { + Scalar minPvValue = eclGrid->getMinpvValue(); + for (int aboveElemCartIdx = static_cast(cartElemIdx) - nx*ny; + aboveElemCartIdx >= 0; + aboveElemCartIdx -= nx*ny) + { + if (porvData[aboveElemCartIdx] >= minPvValue) + // the cartesian element above exhibits a pore volume which larger or + // equal to the minimum one + break; - for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) - porosity_[dofIdx] *= ntgData[gridManager.cartesianIndex(dofIdx)]; - } + Scalar aboveElemVolume = eclGrid->getCellVolume(aboveElemCartIdx); + if (actnumData[aboveElemCartIdx] == 0 && aboveElemVolume > 1e-3) + // stop at explicitly disabled elements, but only if their volume is + // greater than 10^-3 m^3 + break; - // apply the MULTPV keyword to the porosity - if (props.hasDeckDoubleGridProperty("MULTPV")) { - const std::vector &multpvData = - props.getDoubleGridProperty("MULTPV").getData(); - - for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) - porosity_[dofIdx] *= multpvData[gridManager.cartesianIndex(dofIdx)]; - } - - // overwrite the porosity using the PORV keyword for the elements for which PORV - // is defined... - if (props.hasDeckDoubleGridProperty("PORV")) { - const std::vector &porvData = - props.getDoubleGridProperty("PORV").getData(); - - for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) { - unsigned cartesianElemIdx = gridManager.cartesianIndex(dofIdx); - if (std::isfinite(porvData[cartesianElemIdx])) { - Scalar dofVolume = this->simulator().model().dofTotalVolume(dofIdx); - porosity_[dofIdx] = porvData[cartesianElemIdx]/dofVolume; + poreVolume += porvData[aboveElemCartIdx]; } } - } - // the fluid-matrix interactions for ECL problems are dealt with by a separate class + // we define the porosity as the accumulated pore volume divided by the + // geometric volume of the element. Note that -- in pathetic cases -- it can + // be larger than 1.0! + Scalar dofVolume = this->simulator().model().dofTotalVolume(dofIdx); + assert(dofVolume > 0.0); + porosity_[dofIdx] = poreVolume/dofVolume; + } + //////////////////////////////// + + //////////////////////////////// + // fluid-matrix interactions (saturation functions; relperm/capillary pressure) std::vector compressedToCartesianElemIdx(numDof); for (unsigned elemIdx = 0; elemIdx < numDof; ++elemIdx) compressedToCartesianElemIdx[elemIdx] = gridManager.cartesianIndex(elemIdx); materialLawManager_ = std::make_shared(); materialLawManager_->initFromDeck(deck, eclState, compressedToCartesianElemIdx); + //////////////////////////////// } void initFluidSystem_()