diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 8fbc7206c..845ae35f1 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -113,6 +113,7 @@ class EclOutputBlackOilModule : public EclGenericOutputBlackoilModule; using MaterialLaw = GetPropType; using MaterialLawParams = GetPropType; + using IntensiveQuantities = GetPropType; using FluidSystem = GetPropType; using GridView = GetPropType; using Element = typename GridView::template Codim<0>::Entity; @@ -544,128 +545,10 @@ public: this->viscosity_[gasPhaseIdx][globalDofIdx] = FluidSystem::viscosity(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex()); } - - // Add fluid in Place values - //updateFluidInPlace_(elemCtx, dofIdx); - - // Adding block data - const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx); - for (auto& val : this->blockData_) { - const auto& key = val.first; - assert(key.second > 0); - unsigned int cartesianIdxBlock = key.second - 1; - if (cartesianIdx == cartesianIdxBlock) { - if ((key.first == "BWSAT") || (key.first == "BSWAT")) - val.second = getValue(fs.saturation(waterPhaseIdx)); - else if ((key.first == "BGSAT") || (key.first == "BSGAS")) - val.second = getValue(fs.saturation(gasPhaseIdx)); - else if ((key.first == "BOSAT") || (key.first == "BSOIL")) - val.second = getValue(fs.saturation(oilPhaseIdx)); - else if (key.first == "BNSAT") - val.second = intQuants.solventSaturation().value(); - else if ((key.first == "BPR") || (key.first == "BPRESSUR")) { - if (FluidSystem::phaseIsActive(oilPhaseIdx)) - val.second = getValue(fs.pressure(oilPhaseIdx)); - else if (FluidSystem::phaseIsActive(gasPhaseIdx)) - val.second = getValue(fs.pressure(gasPhaseIdx)); - else if (FluidSystem::phaseIsActive(waterPhaseIdx)) - val.second = getValue(fs.pressure(waterPhaseIdx)); - } else if ((key.first == "BTCNFHEA") || (key.first == "BTEMP")) { - if (FluidSystem::phaseIsActive(oilPhaseIdx)) - val.second = getValue(fs.temperature(oilPhaseIdx)); - else if (FluidSystem::phaseIsActive(gasPhaseIdx)) - val.second = getValue(fs.temperature(gasPhaseIdx)); - else if (FluidSystem::phaseIsActive(waterPhaseIdx)) - val.second = getValue(fs.temperature(waterPhaseIdx)); - } else if (key.first == "BWKR" || key.first == "BKRW") - val.second = getValue(intQuants.relativePermeability(waterPhaseIdx)); - else if (key.first == "BGKR" || key.first == "BKRG") - val.second = getValue(intQuants.relativePermeability(gasPhaseIdx)); - else if (key.first == "BOKR" || key.first == "BKRO") - val.second = getValue(intQuants.relativePermeability(oilPhaseIdx)); - else if (key.first == "BKROG") { - const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0); - const auto krog - = MaterialLaw::template relpermOilInOilGasSystem(materialParams, fs); - val.second = getValue(krog); - } else if (key.first == "BKROW") { - const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0); - const auto krow - = MaterialLaw::template relpermOilInOilWaterSystem(materialParams, fs); - val.second = getValue(krow); - } else if (key.first == "BWPC") - val.second = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx)); - else if (key.first == "BGPC") - val.second = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx)); - else if (key.first == "BWPR") - val.second = getValue(fs.pressure(waterPhaseIdx)); - else if (key.first == "BGPR") - val.second = getValue(fs.pressure(gasPhaseIdx)); - else if (key.first == "BVWAT" || key.first == "BWVIS") - val.second = getValue(fs.viscosity(waterPhaseIdx)); - else if (key.first == "BVGAS" || key.first == "BGVIS") - val.second = getValue(fs.viscosity(gasPhaseIdx)); - else if (key.first == "BVOIL" || key.first == "BOVIS") - val.second = getValue(fs.viscosity(oilPhaseIdx)); - else if ((key.first == "BRPV") || (key.first == "BOPV") || (key.first == "BWPV") - || (key.first == "BGPV")) { - if (key.first == "BRPV") { - val.second = 1.0; - } else if (key.first == "BOPV") { - val.second = getValue(fs.saturation(oilPhaseIdx)); - } else if (key.first == "BWPV") { - val.second = getValue(fs.saturation(waterPhaseIdx)); - } else { - val.second = getValue(fs.saturation(gasPhaseIdx)); - } - - // Include active pore-volume. - val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx) - * getValue(intQuants.porosity()); - } else if (key.first == "BRS") - val.second = getValue(fs.Rs()); - else if (key.first == "BRV") - val.second = getValue(fs.Rv()); - else if ((key.first == "BOIP") || (key.first == "BOIPL") || (key.first == "BOIPG") - || (key.first == "BGIP") || (key.first == "BGIPL") || (key.first == "BGIPG") - || (key.first == "BWIP")) { - if ((key.first == "BOIP") || (key.first == "BOIPL")) { - val.second = getValue(fs.invB(oilPhaseIdx)) * getValue(fs.saturation(oilPhaseIdx)); - - if (key.first == "BOIP") { - val.second += getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx)) - * getValue(fs.saturation(gasPhaseIdx)); - } - } else if (key.first == "BOIPG") { - val.second = getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx)) - * getValue(fs.saturation(gasPhaseIdx)); - } else if ((key.first == "BGIP") || (key.first == "BGIPG")) { - val.second = getValue(fs.invB(gasPhaseIdx)) * getValue(fs.saturation(gasPhaseIdx)); - - if (key.first == "BGIP") { - val.second += getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx)) - * getValue(fs.saturation(oilPhaseIdx)); - } - } else if (key.first == "BGIPL") { - val.second = getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx)) - * getValue(fs.saturation(oilPhaseIdx)); - } else { // BWIP - val.second = getValue(fs.invB(waterPhaseIdx)) * getValue(fs.saturation(waterPhaseIdx)); - } - - // Include active pore-volume. - val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx) - * getValue(intQuants.porosity()); - } else { - std::string logstring = "Keyword '"; - logstring.append(key.first); - logstring.append("' is unhandled for output to file."); - OpmLog::warning("Unhandled output keyword", logstring); - } - } - } + // Adding Well RFT data + const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx); if (this->oilConnectionPressures_.count(cartesianIdx) > 0) { this->oilConnectionPressures_[cartesianIdx] = getValue(fs.pressure(oilPhaseIdx)); } @@ -696,9 +579,20 @@ public: = tracerModel.tracerConcentration(tracerIdx, globalDofIdx); } } + } + } - // flows - if (!problem.model().linearizer().getFlowsInfo().empty()) { + void processElementFlows(const ElementContext& elemCtx) + { + OPM_TIMEBLOCK_LOCAL(processElementBlockData); + if (!std::is_same>::value) + return; + + const auto& problem = elemCtx.simulator().problem(); + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { + + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + if (!problem.model().linearizer().getFlowsInfo().empty()) { const auto& flowsInf = problem.model().linearizer().getFlowsInfo(); auto flowsInfos = flowsInf[globalDofIdx]; for (auto& flowsInfo : flowsInfos) { @@ -833,6 +727,137 @@ public: } } + void processElementBlockData(const ElementContext& elemCtx) + { + OPM_TIMEBLOCK_LOCAL(processElementBlockData); + if (!std::is_same>::value) + return; + + const auto& problem = elemCtx.simulator().problem(); + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { + // Adding block data + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx); + const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + for (auto& val : this->blockData_) { + const auto& key = val.first; + assert(key.second > 0); + unsigned int cartesianIdxBlock = key.second - 1; + if (cartesianIdx == cartesianIdxBlock) { + if ((key.first == "BWSAT") || (key.first == "BSWAT")) + val.second = getValue(fs.saturation(waterPhaseIdx)); + else if ((key.first == "BGSAT") || (key.first == "BSGAS")) + val.second = getValue(fs.saturation(gasPhaseIdx)); + else if ((key.first == "BOSAT") || (key.first == "BSOIL")) + val.second = getValue(fs.saturation(oilPhaseIdx)); + else if (key.first == "BNSAT") + val.second = intQuants.solventSaturation().value(); + else if ((key.first == "BPR") || (key.first == "BPRESSUR")) { + if (FluidSystem::phaseIsActive(oilPhaseIdx)) + val.second = getValue(fs.pressure(oilPhaseIdx)); + else if (FluidSystem::phaseIsActive(gasPhaseIdx)) + val.second = getValue(fs.pressure(gasPhaseIdx)); + else if (FluidSystem::phaseIsActive(waterPhaseIdx)) + val.second = getValue(fs.pressure(waterPhaseIdx)); + } else if ((key.first == "BTCNFHEA") || (key.first == "BTEMP")) { + if (FluidSystem::phaseIsActive(oilPhaseIdx)) + val.second = getValue(fs.temperature(oilPhaseIdx)); + else if (FluidSystem::phaseIsActive(gasPhaseIdx)) + val.second = getValue(fs.temperature(gasPhaseIdx)); + else if (FluidSystem::phaseIsActive(waterPhaseIdx)) + val.second = getValue(fs.temperature(waterPhaseIdx)); + } else if (key.first == "BWKR" || key.first == "BKRW") + val.second = getValue(intQuants.relativePermeability(waterPhaseIdx)); + else if (key.first == "BGKR" || key.first == "BKRG") + val.second = getValue(intQuants.relativePermeability(gasPhaseIdx)); + else if (key.first == "BOKR" || key.first == "BKRO") + val.second = getValue(intQuants.relativePermeability(oilPhaseIdx)); + else if (key.first == "BKROG") { + const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0); + const auto krog + = MaterialLaw::template relpermOilInOilGasSystem(materialParams, fs); + val.second = getValue(krog); + } else if (key.first == "BKROW") { + const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0); + const auto krow + = MaterialLaw::template relpermOilInOilWaterSystem(materialParams, fs); + val.second = getValue(krow); + } else if (key.first == "BWPC") + val.second = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx)); + else if (key.first == "BGPC") + val.second = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx)); + else if (key.first == "BWPR") + val.second = getValue(fs.pressure(waterPhaseIdx)); + else if (key.first == "BGPR") + val.second = getValue(fs.pressure(gasPhaseIdx)); + else if (key.first == "BVWAT" || key.first == "BWVIS") + val.second = getValue(fs.viscosity(waterPhaseIdx)); + else if (key.first == "BVGAS" || key.first == "BGVIS") + val.second = getValue(fs.viscosity(gasPhaseIdx)); + else if (key.first == "BVOIL" || key.first == "BOVIS") + val.second = getValue(fs.viscosity(oilPhaseIdx)); + else if ((key.first == "BRPV") || (key.first == "BOPV") || (key.first == "BWPV") + || (key.first == "BGPV")) { + if (key.first == "BRPV") { + val.second = 1.0; + } else if (key.first == "BOPV") { + val.second = getValue(fs.saturation(oilPhaseIdx)); + } else if (key.first == "BWPV") { + val.second = getValue(fs.saturation(waterPhaseIdx)); + } else { + val.second = getValue(fs.saturation(gasPhaseIdx)); + } + + // Include active pore-volume. + val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx) + * getValue(intQuants.porosity()); + } else if (key.first == "BRS") + val.second = getValue(fs.Rs()); + else if (key.first == "BRV") + val.second = getValue(fs.Rv()); + else if ((key.first == "BOIP") || (key.first == "BOIPL") || (key.first == "BOIPG") + || (key.first == "BGIP") || (key.first == "BGIPL") || (key.first == "BGIPG") + || (key.first == "BWIP")) { + if ((key.first == "BOIP") || (key.first == "BOIPL")) { + val.second = getValue(fs.invB(oilPhaseIdx)) * getValue(fs.saturation(oilPhaseIdx)); + + if (key.first == "BOIP") { + val.second += getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx)) + * getValue(fs.saturation(gasPhaseIdx)); + } + } else if (key.first == "BOIPG") { + val.second = getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx)) + * getValue(fs.saturation(gasPhaseIdx)); + } else if ((key.first == "BGIP") || (key.first == "BGIPG")) { + val.second = getValue(fs.invB(gasPhaseIdx)) * getValue(fs.saturation(gasPhaseIdx)); + + if (key.first == "BGIP") { + val.second += getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx)) + * getValue(fs.saturation(oilPhaseIdx)); + } + } else if (key.first == "BGIPL") { + val.second = getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx)) + * getValue(fs.saturation(oilPhaseIdx)); + } else { // BWIP + val.second = getValue(fs.invB(waterPhaseIdx)) * getValue(fs.saturation(waterPhaseIdx)); + } + + // Include active pore-volume. + val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx) + * getValue(intQuants.porosity()); + } else { + std::string logstring = "Keyword '"; + logstring.append(key.first); + logstring.append("' is unhandled for output to file."); + OpmLog::warning("Unhandled output keyword", logstring); + } + } + } + } + } + + /*! * \brief Capture connection fluxes, particularly to account for inter-region flows. * @@ -983,6 +1008,9 @@ public: updateFluidInPlace_(elemCtx, dofIdx); } } + void updateFluidInPlace(unsigned globalDofIdx,const IntensiveQuantities& intQuants, double totVolume){ + this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume); + } private: bool isDefunctParallelWell(std::string wname) const override { @@ -994,12 +1022,19 @@ private: return candidate == parallelWells.end() || *candidate != value; } - void updateFluidInPlace_(const ElementContext& elemCtx, unsigned dofIdx) + void updateFluidInPlace_(const ElementContext& elemCtx, unsigned dofIdx){ + const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0); + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx); + this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume); + } + + void updateFluidInPlace_(unsigned globalDofIdx,const IntensiveQuantities& intQuants, double totVolume) { OPM_TIMEBLOCK_LOCAL(updateFluidInPlace); - const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0); + //const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0); const auto& fs = intQuants.fluidState(); - unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + //unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); // Fluid in Place calculations @@ -1010,7 +1045,7 @@ private: // PORV, MINPV and friends). Also note that because of this, the porosity // returned by the intensive quantities can be outside of the physical // range [0, 1] in pathetic cases. - const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx); + //const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx); const double pv = totVolume * intQuants.porosity().value(); if (!this->pressureTimesHydrocarbonVolume_.empty() && !this->pressureTimesPoreVolume_.empty()) { @@ -1260,7 +1295,7 @@ private: } const Simulator& simulator_; -}; + }; } // namespace Opm diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index 0e5cc9c7e..120124cf9 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -541,13 +541,31 @@ private: eclOutputModule_->processElement(elemCtx); } } + if(!simulator_.model().linearizer().getFlowsInfo().empty()){ + OPM_TIMEBLOCK(prepareFlowsData); + for (const auto& elem : elements(gridView)) { + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + eclOutputModule_->processElementFlows(elemCtx); + } + } { - OPM_TIMEBLOCK(prepareFluidInPlace); + OPM_TIMEBLOCK(prepareBlockData); for (const auto& elem : elements(gridView)) { elemCtx.updatePrimaryStencil(elem); elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - - eclOutputModule_->updateFluidInPlace(elemCtx); + eclOutputModule_->processElementBlockData(elemCtx); + } + } + { + OPM_TIMEBLOCK(prepareFluidInPlace); +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int dofIdx=0; dofIdx < numElements; ++dofIdx){ + const auto& intQuants = *(simulator_.model().cachedIntensiveQuantities(dofIdx, /*timeIdx=*/0)); + const auto totVolume = simulator_.model().dofTotalVolume(dofIdx); + eclOutputModule_->updateFluidInPlace(dofIdx, intQuants, totVolume); } } eclOutputModule_->validateLocalData();