From ecec83349b85ff090c5ce4ed7218ed58095ed008 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 14 Feb 2022 12:30:32 +0100 Subject: [PATCH] Capture Bulk Connection Fluxes to Accumlate Inter-Region Flows This commit adds a new grid traversal that computes fluxes-presently surface level component fluxes-for all bulk connections on the current MPI rank. We aggregate those fluxes, if applicable, into a container for inter-region flows, but this support could be extended to capturing the full 3D vector flow rates for restart output if needed. --- ebos/eclgenericoutputblackoilmodule.cc | 27 +++- ebos/eclgenericoutputblackoilmodule.hh | 5 + ebos/ecloutputblackoilmodule.hh | 190 +++++++++++++++++++++++++ ebos/eclwriter.hh | 44 +++++- 4 files changed, 263 insertions(+), 3 deletions(-) diff --git a/ebos/eclgenericoutputblackoilmodule.cc b/ebos/eclgenericoutputblackoilmodule.cc index 9765a89d9..0aec42c72 100644 --- a/ebos/eclgenericoutputblackoilmodule.cc +++ b/ebos/eclgenericoutputblackoilmodule.cc @@ -18,6 +18,7 @@ */ #include + #include #include @@ -33,11 +34,16 @@ #include #include -#include +#include +#include #include #include #include #include +#include +#include + +#include namespace { @@ -59,6 +65,24 @@ std::string EclString(Opm::Inplace::Phase phase) { } } + std::size_t numCells(const Opm::EclipseState& eclState) + { + return eclState.fieldProps().get_int("FIPNUM").size(); + } + + std::vector + defineInterRegionFlowArrays(const Opm::EclipseState& eclState, + const Opm::SummaryConfig& summaryConfig) + { + auto regions = std::vector{}; + + const auto& fprops = eclState.fieldProps(); + for (const auto& arrayName : summaryConfig.fip_regions_interreg_flow()) { + regions.push_back({ arrayName, std::cref(fprops.get_int(arrayName)) }); + } + + return regions; + } } namespace Opm { @@ -82,6 +106,7 @@ EclGenericOutputBlackoilModule(const EclipseState& eclState, , schedule_(schedule) , summaryConfig_(summaryConfig) , summaryState_(summaryState) + , interRegionFlows_(numCells(eclState), defineInterRegionFlowArrays(eclState, summaryConfig)) , enableEnergy_(enableEnergy) , enableTemperature_(enableTemperature) , enableSolvent_(enableSolvent) diff --git a/ebos/eclgenericoutputblackoilmodule.hh b/ebos/eclgenericoutputblackoilmodule.hh index 3c17c3498..915a24885 100644 --- a/ebos/eclgenericoutputblackoilmodule.hh +++ b/ebos/eclgenericoutputblackoilmodule.hh @@ -43,6 +43,8 @@ #endif #include +#include + namespace Opm { namespace data { class Solution; } @@ -370,6 +372,9 @@ protected: const Schedule& schedule_; const SummaryConfig& summaryConfig_; const SummaryState& summaryState_; + + EclInterRegFlowMap interRegionFlows_; + bool enableEnergy_; bool enableTemperature_; diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index b3fe26095..0306ce9f3 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -31,6 +31,8 @@ #include #include #include +#include +#include #include @@ -576,6 +578,94 @@ public: } } + /*! + * \brief Capture connection fluxes, particularly to account for inter-region flows. + * + * \tparam ActiveIndex Callable type, typically a lambda, that enables + * retrieving the active index, on the local MPI rank, of a + * particular cell/element. Must support a function call operator of + * the form + \code + int operator()(const Element& elem) const + \endcode + * + * \tparam CartesianIndex Callable type, typically a lambda, that + * enables retrieving the globally unique Cartesian index of a + * particular cell/element given its active index on the local MPI + * rank. Must support a function call operator of the form + \code + int operator()(const int activeIndex) const + \endcode + * + * \param[in] elemCtx Primary lookup structure for per-cell/element + * dynamic information. + * + * \param[in] activeIndex Mapping from cell/elements to linear indices + * on local MPI rank. + * + * \param[in] cartesianIndex Mapping from active index on local MPI rank + * to globally unique Cartesian cell/element index. + */ + template + void processFluxes(const ElementContext& elemCtx, + ActiveIndex&& activeIndex, + CartesianIndex&& cartesianIndex) + { + const auto identifyCell = [&activeIndex, &cartesianIndex](const Element& elem) + -> EclInterRegFlowMap::Cell + { + const auto cellIndex = activeIndex(elem); + + return { + cellIndex, + cartesianIndex(cellIndex), + elem.partitionType() == Dune::InteriorEntity + }; + }; + + const auto timeIdx = 0u; + const auto& stencil = elemCtx.stencil(timeIdx); + const auto numInteriorFaces = elemCtx.numInteriorFaces(timeIdx); + + for (auto scvfIdx = 0*numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) { + const auto& face = stencil.interiorFace(scvfIdx); + const auto left = identifyCell(stencil.element(face.interiorIndex())); + const auto right = identifyCell(stencil.element(face.exteriorIndex())); + + const auto rates = this-> + getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx); + + this->interRegionFlows_.addConnection(left, right, rates); + } + } + + /*! + * \brief Prepare for capturing connection fluxes, particularly to + * account for inter-region flows. + */ + void initializeFluxData() + { + // Inter-region flow rates. Note: ".clear()" prepares to accumulate + // contributions per bulk connection between FIP regions. + this->interRegionFlows_.clear(); + } + + /*! + * \brief Finalize capturing connection fluxes. + */ + void finalizeFluxData() + { + this->interRegionFlows_.compress(); + } + + /*! + * \brief Get read-only access to collection of inter-region flows. + */ + const EclInterRegFlowMap& getInterRegFlows() const + { + return this->interRegionFlows_; + } + template void assignToFluidState(FluidState& fs, unsigned elemIdx) const { @@ -772,6 +862,106 @@ private: } } + /*! + * \brief Compute surface level component flow rates across a single + * intersection. + * + * \param[in] elemCtx Primary lookup structure for per-cell/element + * dynamic information. + * + * \param[in] scvfIdx Linear index of current interior bulk connection. + * + * \param[in] timeIdx Historical time-point at which to evaluate dynamic + * quantities (e.g., reciprocal FVF or dissolved gas concentration). + * Zero for the current time. + * + * \return Surface level component flow rates. + */ + data::InterRegFlowMap::FlowRates + getComponentSurfaceRates(const ElementContext& elemCtx, + const Scalar faceArea, + const std::size_t scvfIdx, + const std::size_t timeIdx) const + { + using Component = data::InterRegFlowMap::Component; + + auto rates = data::InterRegFlowMap::FlowRates{}; + + const auto& extQuant = elemCtx.extensiveQuantities(scvfIdx, timeIdx); + + const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea; + + if (FluidSystem::phaseIsActive(oilPhaseIdx)) { + const auto& up = elemCtx + .intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx); + + using FluidState = std::remove_cv_t>; + + const auto pvtReg = up.pvtRegionIndex(); + + const auto bO = getValue(getInvB_ + (up.fluidState(), oilPhaseIdx, pvtReg)); + + const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx)); + + rates[Component::Oil] += qO; + + if (FluidSystem::phaseIsActive(gasPhaseIdx)) { + const auto Rs = getValue( + BlackOil::getRs_ + (up.fluidState(), pvtReg)); + + rates[Component::Gas] += qO * Rs; + rates[Component::Disgas] += qO * Rs; + } + } + + if (FluidSystem::phaseIsActive(gasPhaseIdx)) { + const auto& up = elemCtx + .intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx); + + using FluidState = std::remove_cv_t>; + + const auto pvtReg = up.pvtRegionIndex(); + + const auto bG = getValue(getInvB_ + (up.fluidState(), gasPhaseIdx, pvtReg)); + + const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx)); + + rates[Component::Gas] += qG; + + if (FluidSystem::phaseIsActive(oilPhaseIdx)) { + const auto Rv = getValue( + BlackOil::getRv_ + (up.fluidState(), pvtReg)); + + rates[Component::Oil] += qG * Rv; + rates[Component::Vapoil] += qG * Rv; + } + } + + if (FluidSystem::phaseIsActive(waterPhaseIdx)) { + const auto& up = elemCtx + .intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx); + + using FluidState = std::remove_cv_t>; + + const auto pvtReg = up.pvtRegionIndex(); + + const auto bW = getValue(getInvB_ + (up.fluidState(), waterPhaseIdx, pvtReg)); + + rates[Component::Water] += + alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx)); + } + + return rates; + } + const Simulator& simulator_; }; diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index 043dd58fa..dee9a66ae 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -38,8 +38,11 @@ #include -#include +#include + +#include #include +#include namespace Opm::Properties { @@ -181,11 +184,12 @@ public: const auto localGroupAndNetworkData = simulator_.problem().wellModel() .groupAndNetworkData(reportStepNum); - const auto localAquiferData = simulator_.problem().aquiferModel().aquiferData(); const auto localWellTestState = simulator_.problem().wellModel().wellTestState(); this->prepareLocalCellData(isSubStep, reportStepNum); + this->captureLocalFluxData(); + if (this->collectToIORank_.isParallel()) this->collectToIORank_.collect({}, eclOutputModule_->getBlockData(), @@ -435,6 +439,42 @@ private: OPM_END_PARALLEL_TRY_CATCH("EclWriter::prepareLocalCellData() failed: ", simulator_.vanguard().grid().comm()) } + void captureLocalFluxData() + { + const auto& gridView = this->simulator_.vanguard().gridView(); + const auto timeIdx = 0u; + + auto elemCtx = ElementContext { this->simulator_ }; + + const auto elemMapper = ElementMapper { gridView, Dune::mcmgElementLayout() }; + const auto activeIndex = [&elemMapper](const Element& e) + { + return elemMapper.index(e); + }; + + const auto cartesianIndex = [this](const int elemIndex) + { + return this->cartMapper_.cartesianIndex(elemIndex); + }; + + this->eclOutputModule_->initializeFluxData(); + + OPM_BEGIN_PARALLEL_TRY_CATCH(); + + for (const auto& elem : elements(gridView, Dune::Partitions::interiorBorder)) { + elemCtx.updateStencil(elem); + elemCtx.updateIntensiveQuantities(timeIdx); + elemCtx.updateExtensiveQuantities(timeIdx); + + this->eclOutputModule_->processFluxes(elemCtx, activeIndex, cartesianIndex); + } + + OPM_END_PARALLEL_TRY_CATCH("EclWriter::captureLocalFluxData() failed: ", + this->simulator_.vanguard().grid().comm()) + + this->eclOutputModule_->finalizeFluxData(); + } + Simulator& simulator_; std::unique_ptr> eclOutputModule_; Scalar restartTimeStepSize_;