From f62498440d9433f1fde9a4c9c6534674a8a6a42b Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 12 Jan 2018 15:32:43 +0100 Subject: [PATCH] Add FIP, INIT and EGRID output Writes INIT and EGRID files initially Adds Fip to summary output if required. Output Fip to log (.PRT) if Opm-logger is set --- ebos/collecttoiorank.hh | 24 +- ebos/ecloutputblackoilmodule.hh | 375 ++++++++++++++++++++++++++++++-- ebos/eclproblem.hh | 12 +- ebos/eclwriter.hh | 182 +++++++++++++++- 4 files changed, 551 insertions(+), 42 deletions(-) diff --git a/ebos/collecttoiorank.hh b/ebos/collecttoiorank.hh index 6d4f94722..9c14eb663 100644 --- a/ebos/collecttoiorank.hh +++ b/ebos/collecttoiorank.hh @@ -93,16 +93,19 @@ namespace Ewoms IndexMapType& localIndexMap_; IndexMapStorageType& indexMaps_; std::map< const int, const int > globalPosition_; + std::vector& ranks_; public: DistributeIndexMapping( const std::vector& globalIndex, const std::vector& distributedGlobalIndex, IndexMapType& localIndexMap, - IndexMapStorageType& indexMaps ) + IndexMapStorageType& indexMaps, + std::vector& ranks) : distributedGlobalIndex_( distributedGlobalIndex ), localIndexMap_( localIndexMap ), indexMaps_( indexMaps ), - globalPosition_() + globalPosition_(), + ranks_(ranks) { const size_t size = globalIndex.size(); // create mapping globalIndex --> localIndex @@ -114,6 +117,7 @@ namespace Ewoms // we need to create a mapping from local to global if( ! indexMaps_.empty() ) { + ranks_.resize( size, -1); // for the ioRank create a localIndex to index in global state map IndexMapType& indexMap = indexMaps_.back(); const size_t localSize = localIndexMap_.size(); @@ -122,6 +126,7 @@ namespace Ewoms { const int id = distributedGlobalIndex_[ localIndexMap_[ i ] ]; indexMap[ i ] = globalPosition_[ id ] ; + ranks_[ indexMap[ i ] ] = ioRank; } } } @@ -160,6 +165,7 @@ namespace Ewoms buffer.read( globalId ); assert( globalPosition_.find( globalId ) != globalPosition_.end() ); indexMap[ index ] = globalPosition_[ globalId ]; + ranks_[ indexMap[ index ] ] = link + 1; } } }; @@ -265,7 +271,7 @@ namespace Ewoms indexMaps_.resize( comm.size() ); // distribute global id's to io rank for later association of dof's - DistributeIndexMapping distIndexMapping( globalCartesianIndex_, distributedCartesianIndex, localIndexMap_, indexMaps_ ); + DistributeIndexMapping distIndexMapping( globalCartesianIndex_, distributedCartesianIndex, localIndexMap_, indexMaps_, globalRanks_); toIORankComm_.exchange( distIndexMapping ); } } @@ -435,14 +441,14 @@ namespace Ewoms return toIORankComm_.size() > 1; } - int localIdxToGlobalIdx(const unsigned localIdx) { - + int localIdxToGlobalIdx(const unsigned localIdx) const + { if ( ! isParallel() ) { return localIdx; } // the last indexMap is the local one - IndexMapType& indexMap = indexMaps_.back(); + const IndexMapType& indexMap = indexMaps_.back(); if( indexMap.empty() ) OPM_THROW(std::logic_error,"index map is not created on this rank"); @@ -454,11 +460,17 @@ namespace Ewoms size_t numCells () const { return globalCartesianIndex_.size(); } + const std::vector& globalRanks() const + { + return globalRanks_; + } + protected: P2PCommunicatorType toIORankComm_; IndexMapType globalCartesianIndex_; IndexMapType localIndexMap_; IndexMapStorageType indexMaps_; + std::vector globalRanks_; Opm::data::Solution globalCellData_; }; diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 9e6653038..1bfa56f82 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -34,6 +34,7 @@ #include #include +#include #include #include @@ -42,6 +43,8 @@ #include +#define ENUM_TO_STR(ENUM) std::string(#ENUM) + namespace Ewoms { namespace Properties { // create new type tag for the Ecl-output @@ -81,10 +84,28 @@ class EclOutputBlackOilModule typedef std::vector ScalarBuffer; + struct FIPDataType { + + enum FipId { + WIP = waterPhaseIdx, + OIP = oilPhaseIdx, + GIP = gasPhaseIdx, + OIPL = 3, + OIPG = 4, + GIPL = 5, + GIPG = 6, + PV = 7, //< Pore volume + PAV = 8 + }; + static const int fipValues = PAV + 1 ; + }; + public: EclOutputBlackOilModule(const Simulator& simulator) : simulator_(simulator) - { } + { + createLocalFipnum_(); + } /*! * \brief Allocate memory for the scalar fields we would like to @@ -101,37 +122,61 @@ public: keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum); } + const Opm::SummaryConfig summaryConfig = simulator_.gridManager().summaryConfig(); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) continue; - saturation_[phaseIdx].resize(bufferSize,0.0); + if (!substep || (phaseIdx == waterPhaseIdx && summaryConfig.require3DField("SWAT") ) + || (phaseIdx == gasPhaseIdx && summaryConfig.require3DField("SGAS") ) ) + saturation_[phaseIdx].resize(bufferSize,0.0); } - oilPressure_.resize(bufferSize,0.0); - temperature_.resize(bufferSize,0.0); + if (!substep || summaryConfig.require3DField("PRESSURE")) + oilPressure_.resize(bufferSize,0.0); + + if (!substep || summaryConfig.require3DField("TEMP")) + temperature_.resize(bufferSize,0.0); // Output the same as legacy // TODO: Only needed if DISGAS or VAPOIL - if (true) - rs_.resize(bufferSize,0.0); - if (true) - rv_.resize(bufferSize,0.0); + if (true) { + if (!substep || summaryConfig.require3DField("RS")) + rs_.resize(bufferSize,0.0); + } + if (true) { + if (!substep || summaryConfig.require3DField("RV")) + rv_.resize(bufferSize,0.0); + } if (GET_PROP_VALUE(TypeTag, EnableSolvent)) { - sSol_.resize(bufferSize,0.0); + if (!substep || summaryConfig.require3DField("SSOL")) + sSol_.resize(bufferSize,0.0); } if (GET_PROP_VALUE(TypeTag, EnablePolymer)) { - cPolymer_.resize(bufferSize,0.0); + if (!substep || summaryConfig.require3DField("POLYMER")) + cPolymer_.resize(bufferSize,0.0); } + // Fluid in place + for (int i = 0; i::type>::type FluidState; unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); unsigned pvtRegionIdx = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex(); @@ -438,6 +484,62 @@ public: intQuants.pvtRegionIndex()); } + // FIP + + // calculate the pore volume of the current cell. Note that the porosity + // returned by the intensive quantities is defined as the ratio of pore + // space to total cell volume and includes all pressure dependent (-> + // rock compressibility) and static modifiers (MULTPV, MULTREGP, NTG, + // 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 double pv = + elemCtx.simulator().model().dofTotalVolume(globalDofIdx) + * intQuants.porosity().value(); + if (fip_[FIPDataType::PV].size() > 0) + fip_[FIPDataType::PV][globalDofIdx] = pv; + + Scalar fip[FluidSystem::numPhases]; + for (unsigned phase = 0; phase < FluidSystem::numPhases; ++phase) { + if (!FluidSystem::phaseIsActive(phase)) { + continue; + } + + const double b = Toolbox::value(fs.invB(phase)); + const double s = Toolbox::value(fs.saturation(phase)); + fip[phase] = b * s * pv; + } + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && fip_[FIPDataType::OIP].size() > 0) + fip_[FIPDataType::OIP][globalDofIdx] = fip[FluidSystem::oilPhaseIdx]; + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && fip_[FIPDataType::GIP].size() > 0) + fip_[FIPDataType::GIP][globalDofIdx] = fip[FluidSystem::gasPhaseIdx]; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && fip_[FIPDataType::WIP].size() > 0) + fip_[FIPDataType::WIP][globalDofIdx] = fip[FluidSystem::waterPhaseIdx]; + + // Store the pure oil and gas FIP + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && fip_[FIPDataType::OIPL].size() > 0) + fip_[FIPDataType::OIPL][globalDofIdx] = fip[FluidSystem::oilPhaseIdx]; + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && fip_[FIPDataType::GIPG].size() > 0) + fip_[FIPDataType::GIPG][globalDofIdx] = fip[FluidSystem::gasPhaseIdx]; + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + // Gas dissolved in oil and vaporized oil + Scalar gipl = Toolbox::value(fs.Rs()) * fip[FluidSystem::oilPhaseIdx]; + Scalar oipg = Toolbox::value(fs.Rv()) * fip[FluidSystem::gasPhaseIdx]; + if (fip_[FIPDataType::GIPG].size() > 0) + fip_[FIPDataType::GIPL][globalDofIdx] = gipl; + if (fip_[FIPDataType::OIPG].size() > 0) + fip_[FIPDataType::OIPG][globalDofIdx] = oipg; + + // Add dissolved gas and vaporized oil to total FIP + if (fip_[FIPDataType::OIP].size() > 0) + fip_[FIPDataType::OIP][globalDofIdx] += oipg; + if (fip_[FIPDataType::GIP].size() > 0) + fip_[FIPDataType::GIP][globalDofIdx] += gipl; + } + } } @@ -519,7 +621,7 @@ public: return; if ( oilPressure_.size() > 0 ) { - sol.insert( "PRESSURE", Opm::UnitSystem::measure::pressure, std::move(oilPressure_), Opm::data::TargetType::RESTART_SOLUTION); + sol.insert( "PRESSURE", Opm::UnitSystem::measure::pressure, oilPressure_, Opm::data::TargetType::RESTART_SOLUTION); } if ( temperature_.size() > 0 ) { @@ -527,10 +629,10 @@ public: } if( FluidSystem::phaseIsActive(waterPhaseIdx) && saturation_[waterPhaseIdx].size() > 0 ) { - sol.insert( "SWAT", Opm::UnitSystem::measure::identity, std::move(saturation_[waterPhaseIdx]), Opm::data::TargetType::RESTART_SOLUTION ); + sol.insert( "SWAT", Opm::UnitSystem::measure::identity, saturation_[waterPhaseIdx], Opm::data::TargetType::RESTART_SOLUTION ); } if( FluidSystem::phaseIsActive(gasPhaseIdx) && saturation_[gasPhaseIdx].size() > 0) { - sol.insert( "SGAS", Opm::UnitSystem::measure::identity, std::move(saturation_[gasPhaseIdx]), Opm::data::TargetType::RESTART_SOLUTION ); + sol.insert( "SGAS", Opm::UnitSystem::measure::identity, saturation_[gasPhaseIdx], Opm::data::TargetType::RESTART_SOLUTION ); } if ( gasDissolutionFactor_.size() > 0 ) { @@ -614,7 +716,79 @@ public: if (bubblePointPressure_.size() > 0) sol.insert ("PBUB", Opm::UnitSystem::measure::pressure, std::move(bubblePointPressure_), Opm::data::TargetType::RESTART_AUXILIARY); -} + // Summary FIP output + // Fluid in place + for (int i = 0; i 0) { + sol.insert(stringOfEnumIndex_(i), + Opm::UnitSystem::measure::volume, + fip_[i] , + Opm::data::TargetType::SUMMARY); + } + } + + } + + // write Fip to output log + void outputFIPLog() { + + const auto& comm = simulator_.gridView().comm(); + size_t ntFip = *std::max_element(fipnum_.begin(), fipnum_.end()); + ntFip = comm.max(ntFip); + + // sum values over each region + ScalarBuffer regionValues[FIPDataType::fipValues]; + for (int i = 0; i fieldNum(ntFip, 1); + ScalarBuffer totalValues(FIPDataType::fipValues,0.0); + for (int i = 0; i& fipnum_global = simulator_.gridManager().eclState().get3DProperties().getIntGridProperty("FIPNUM").getData(); + // Get compressed cell fipnum. + const auto& gridView = simulator_.gridManager().gridView(); + unsigned numElements = gridView.size(/*codim=*/0); + fipnum_.resize(numElements); + if (fipnum_global.empty()) { + std::fill(fipnum_.begin(), fipnum_.end(), 0); + } else { + for (size_t elemIdx = 0; elemIdx < numElements; ++elemIdx) { + fipnum_[elemIdx] = fipnum_global[simulator_.gridManager().cartesianIndex( elemIdx )]; + } + } + } + + // Sum Fip values over regions. + ScalarBuffer FIPTotals_(const ScalarBuffer& fip, std::vector& regionId, size_t maxNumberOfRegions, bool commSum = true) + { + ScalarBuffer totals(maxNumberOfRegions, 0.0); + assert(regionId.size() == fip.size()); + for (size_t j = 0; j < regionId.size(); ++j) { + const int regionIdx = regionId[j] - 1; + // the cell is not attributed to any region. ignore it! + if (regionIdx < 0) + continue; + + assert(regionIdx < static_cast(maxNumberOfRegions)); + totals[regionIdx] += fip[j]; + } + if (commSum) { + const auto& comm = simulator_.gridView().comm(); + for (size_t i = 0; i < maxNumberOfRegions; ++i) + totals[i] = comm.sum(totals[i]); + } + + return totals; + } + + // computes the hydrocarbon volume weighted averaged pressure + // of a region (reg) + // if reg == -1, the field average is computed. + Scalar FipPav_(int reg) + { + Scalar pPvSum = 0.0; + Scalar pvSum = 0.0; + Scalar pPvHydrocarbonSum = 0.0; + Scalar pvHydrocarbonSum = 0.0; + size_t numElem = oilPressure_.size(); + for (size_t elem = 0; elem < numElem; ++elem) { + if(static_cast(fipnum_[elem]) == reg || -1 == reg) + { + Scalar hydrocarbon = 0.0; + if (FluidSystem::phaseIsActive(oilPhaseIdx)) + hydrocarbon += saturation_[oilPhaseIdx][elem]; + if (FluidSystem::phaseIsActive(gasPhaseIdx)) + hydrocarbon += saturation_[gasPhaseIdx][elem]; + + pPvSum += oilPressure_[elem] * fip_[FIPDataType::PV][elem]; + pvSum += fip_[FIPDataType::PV][elem]; + pPvHydrocarbonSum += pPvSum * hydrocarbon; + pvHydrocarbonSum += pvSum * hydrocarbon; + } + } + const auto& comm = simulator_.gridView().comm(); + pPvSum = comm.sum(pPvSum); + pvSum = comm.sum(pvSum); + pPvHydrocarbonSum = comm.sum(pPvHydrocarbonSum); + pvHydrocarbonSum = comm.sum(pvHydrocarbonSum); + + if (pvHydrocarbonSum > 1e-10) + return pPvHydrocarbonSum / pvHydrocarbonSum; + + // return the porevolume weighted pressure if no hydrocarbon + return pPvSum / pvSum;; + } + + void FIPUnitConvert_(const Opm::UnitSystem& units, + ScalarBuffer& fip) + { + if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) { + fip[FIPDataType::WIP] = Opm::unit::convert::to(fip[FIPDataType::WIP], Opm::unit::stb); + fip[FIPDataType::OIP] = Opm::unit::convert::to(fip[FIPDataType::OIP], Opm::unit::stb); + fip[FIPDataType::OIPL] = Opm::unit::convert::to(fip[FIPDataType::OIPL], Opm::unit::stb); + fip[FIPDataType::OIPG] = Opm::unit::convert::to(fip[FIPDataType::OIPG], Opm::unit::stb); + fip[FIPDataType::GIP] = Opm::unit::convert::to(fip[FIPDataType::GIP], 1000*Opm::unit::cubic(Opm::unit::feet)); + fip[FIPDataType::GIPL] = Opm::unit::convert::to(fip[FIPDataType::GIPL], 1000*Opm::unit::cubic(Opm::unit::feet)); + fip[FIPDataType::GIPG] = Opm::unit::convert::to(fip[FIPDataType::GIPG], 1000*Opm::unit::cubic(Opm::unit::feet)); + fip[FIPDataType::PV] = Opm::unit::convert::to(fip[FIPDataType::PV], Opm::unit::stb); + fip[FIPDataType::PAV] = Opm::unit::convert::to(fip[FIPDataType::PAV], Opm::unit::psia); + } + else if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_METRIC) { + fip[FIPDataType::PAV] = Opm::unit::convert::to(fip[FIPDataType::PAV], Opm::unit::barsa); + } + else { + OPM_THROW(std::runtime_error, "Unsupported unit type for fluid in place output."); + } + } + + void outputRegionFluidInPlace_(const ScalarBuffer& oip, const ScalarBuffer& cip, const Opm::UnitSystem& units, const int reg) + { + std::ostringstream ss; + if (!reg) { + ss << " ===================================================\n" + << " : Field Totals :\n"; + } else { + ss << " ===================================================\n" + << " : FIPNUM report region " + << std::setw(2) << reg << " :\n"; + } + if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_METRIC) { + ss << " : PAV =" << std::setw(14) << cip[FIPDataType::PAV] << " BARSA :\n" + << std::fixed << std::setprecision(0) + << " : PORV =" << std::setw(14) << cip[FIPDataType::PV] << " RM3 :\n"; + if (!reg) { + ss << " : Pressure is weighted by hydrocarbon pore volume :\n" + << " : Porv volumes are taken at reference conditions :\n"; + } + ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n"; + } + if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) { + ss << " : PAV =" << std::setw(14) << cip[FIPDataType::PAV] << " PSIA :\n" + << std::fixed << std::setprecision(0) + << " : PORV =" << std::setw(14) << cip[FIPDataType::PV] << " RB :\n"; + if (!reg) { + ss << " : Pressure is weighted by hydrocarbon pore volume :\n" + << " : Pore volumes are taken at reference conditions :\n"; + } + ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n"; + } + ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n" + << ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n" + << ":Currently in place :" << std::setw(14) << cip[FIPDataType::OIPL] << std::setw(14) << cip[FIPDataType::OIPG] << std::setw(14) << cip[FIPDataType::OIP] << ":" + << std::setw(13) << cip[FIPDataType::WIP] << " :" << std::setw(14) << (cip[FIPDataType::GIPG]) << std::setw(14) << cip[FIPDataType::GIPL] << std::setw(14) << cip[FIPDataType::GIP] << ":\n" + << ":------------------------:------------------------------------------:----------------:------------------------------------------:\n" + << ":Originally in place :" << std::setw(14) << oip[FIPDataType::OIPL] << std::setw(14) << oip[FIPDataType::OIPG] << std::setw(14) << oip[FIPDataType::OIP] << ":" + << std::setw(13) << oip[FIPDataType::WIP] << " :" << std::setw(14) << oip[FIPDataType::GIPG] << std::setw(14) << oip[FIPDataType::GIPL] << std::setw(14) << oip[FIPDataType::GIP] << ":\n" + << ":========================:==========================================:================:==========================================:\n"; + Opm::OpmLog::note(ss.str()); + } + + std::string stringOfEnumIndex_(int i) { + typedef typename FIPDataType::FipId FipId; + switch( static_cast(i) ) + { + case FIPDataType::WIP: return "WIP"; break; + case FIPDataType::OIP: return "OIP"; break; + case FIPDataType::GIP: return "GIP"; break; + case FIPDataType::OIPL:return "OIPL"; break; + case FIPDataType::OIPG:return "OIPG"; break; + case FIPDataType::GIPL:return "GIPL"; break; + case FIPDataType::GIPG:return "GIPG"; break; + case FIPDataType::PV: return "PV"; break; + case FIPDataType::PAV: return "PAV"; break; + } + return "ERROR"; + } + + const Simulator& simulator_; ScalarBuffer saturation_[numPhases]; @@ -794,6 +1127,10 @@ private: ScalarBuffer dewPointPressure_; std::vector failedCellsPb_; std::vector failedCellsPd_; + std::vector fipnum_; + ScalarBuffer fip_[FIPDataType::fipValues]; + ScalarBuffer origTotalValues_; + ScalarBuffer origRegionValues_[FIPDataType::fipValues]; }; } // namespace Ewoms diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index e4323b2f3..1d4bc7ca7 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -460,6 +460,8 @@ public: int numElements = gridView.size(/*codim=*/0); maxPolymerAdsorption_.resize(numElements, 0.0); } + + eclWriter_->writeInit(); } void prefetch(const Element& elem) const @@ -702,11 +704,10 @@ public: } Scalar totalSolverTime = 0.0; Scalar nextstep = this->simulator().timeStepSize(); - Opm::data::Solution fip; - writeOutput(dw, t, false, totalSolverTime, nextstep, fip, verbose); + writeOutput(dw, t, false, totalSolverTime, nextstep, verbose); } - void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep, const Opm::data::Solution& fip, bool verbose = true) + void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep, bool verbose = true) { // use the generic code to prepare the output fields and to // write the desired VTK files. @@ -714,7 +715,7 @@ public: // output using eclWriter if enabled if (eclWriter_) - eclWriter_->writeOutput(dw, t, substep, totalSolverTime, nextstep, fip); + eclWriter_->writeOutput(dw, t, substep, totalSolverTime, nextstep); } /*! @@ -1166,9 +1167,6 @@ public: return initialFluidStates_[globalDofIdx]; } - void setEclIO(std::unique_ptr&& eclIO) - { eclWriter_->setEclIO(std::move(eclIO)); } - const Opm::EclipseIO& eclIO() const { return eclWriter_->eclIO(); } diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index 6c794327c..9d3b61547 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -95,10 +95,10 @@ public: , eclOutputModule_(simulator) , collectToIORank_(simulator_.gridManager()) { - Grid globalGrid = simulator_.gridManager().grid(); - globalGrid.switchToGlobalView(); + globalGrid_ = simulator_.gridManager().grid(); + globalGrid_.switchToGlobalView(); eclIO_.reset(new Opm::EclipseIO(simulator_.gridManager().eclState(), - Opm::UgGridHelpers::createEclipseGrid( globalGrid , simulator_.gridManager().eclState().getInputGrid() ), + Opm::UgGridHelpers::createEclipseGrid( globalGrid_ , simulator_.gridManager().eclState().getInputGrid() ), simulator_.gridManager().schedule(), simulator_.gridManager().summaryConfig())); } @@ -106,16 +106,28 @@ public: ~EclWriter() { } - void setEclIO(std::unique_ptr&& eclIO) - { eclIO_ = std::move(eclIO); } - const Opm::EclipseIO& eclIO() const { return *eclIO_; } + void writeInit() + { +#if !HAVE_OPM_OUTPUT + OPM_THROW(std::runtime_error, + "Opm-output must be available to write ECL output!"); +#else + if (collectToIORank_.isIORank()) { + std::map > integerVectors; + if (collectToIORank_.isParallel()) + integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks()); + eclIO_->writeInitial(computeTrans_(), integerVectors, exportNncStructure_()); + } +#endif + } + /*! * \brief collect and pass data and pass it to eclIO writer */ - void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep, const Opm::data::Solution& fip) + void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep) { #if !HAVE_OPM_OUTPUT OPM_THROW(std::runtime_error, @@ -133,6 +145,9 @@ public: const ElementIterator& elemEndIt = gridView.template end(); for (; elemIt != elemEndIt; ++elemIt) { const Element& elem = *elemIt; + if (elem.partitionType() != Dune::InteriorEntity) + continue; + elemCtx.updatePrimaryStencil(elem); elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); eclOutputModule_.processElement(elemCtx); @@ -140,11 +155,14 @@ public: eclOutputModule_.outputErrorLog(); // collect all data to I/O rank and assign to sol - Opm::data::Solution localCellData = fip; + Opm::data::Solution localCellData; eclOutputModule_.assignToSolution(localCellData); if (collectToIORank_.isParallel()) collectToIORank_.collect(localCellData); + if (!substep) + eclOutputModule_.outputFIPLog(); + // write output on I/O rank if (collectToIORank_.isIORank()) { @@ -170,7 +188,7 @@ public: false); } -#endif + #endif } void restartBegin() @@ -194,7 +212,7 @@ public: unsigned episodeIdx = simulator_.episodeIndex(); const auto& gridView = simulator_.gridManager().gridView(); unsigned numElements = gridView.size(/*codim=*/0); - eclOutputModule_.allocBuffers(numElements, episodeIdx, simulator_.gridManager().eclState().getRestartConfig(), true, false); + eclOutputModule_.allocBuffers(numElements, episodeIdx, simulator_.gridManager().eclState().getRestartConfig(), false, false); auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys); for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { @@ -213,10 +231,154 @@ private: static bool enableEclOutput_() { return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); } + Opm::data::Solution computeTrans_() const + { + const auto& cartMapper = simulator_.gridManager().cartesianIndexMapper(); + const auto& cartDims = cartMapper.cartesianDimensions(); + const int globalSize = cartDims[0]*cartDims[1]*cartDims[2]; + + Opm::data::CellData tranx = {Opm::UnitSystem::measure::transmissibility, std::vector( globalSize ), Opm::data::TargetType::INIT}; + Opm::data::CellData trany = {Opm::UnitSystem::measure::transmissibility, std::vector( globalSize ), Opm::data::TargetType::INIT}; + Opm::data::CellData tranz = {Opm::UnitSystem::measure::transmissibility, std::vector( globalSize ), Opm::data::TargetType::INIT}; + + for (size_t i = 0; i < tranx.data.size(); ++i) { + tranx.data[0] = 0.0; + trany.data[0] = 0.0; + tranz.data[0] = 0.0; + } + + const auto& globalGridView = globalGrid_.leafGridView(); + typedef typename Grid::LeafGridView GridView; + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; + ElementMapper globalElemMapper(globalGridView); + const auto& cartesianCellIdx = globalGrid_.globalCell(); + + const auto* globalTrans = &(simulator_.gridManager().globalTransmissibility()); + if (!collectToIORank_.isParallel()) { + // in the sequential case we must use the transmissibilites defined by + // the problem. (because in the sequential case, the grid manager does + // not compute "global" transmissibilities for performance reasons. in + // the parallel case, the problem's transmissibilities can't be used + // because this object refers to the distributed grid and we need the + // sequential version here.) + globalTrans = &simulator_.problem().eclTransmissibilities(); + } + + auto elemIt = globalGridView.template begin(); + const auto& elemEndIt = globalGridView.template end(); + for (; elemIt != elemEndIt; ++ elemIt) { + const auto& elem = *elemIt; + + auto isIt = globalGridView.ibegin(elem); + const auto& isEndIt = globalGridView.iend(elem); + for (; isIt != isEndIt; ++ isIt) { + const auto& is = *isIt; + + if (!is.neighbor()) + { + continue; // intersection is on the domain boundary + } + + unsigned c1 = globalElemMapper.index(is.inside()); + unsigned c2 = globalElemMapper.index(is.outside()); + + if (c1 > c2) + { + continue; // we only need to handle each connection once, thank you. + } + + + int gc1 = std::min(cartesianCellIdx[c1], cartesianCellIdx[c2]); + int gc2 = std::max(cartesianCellIdx[c1], cartesianCellIdx[c2]); + if (gc2 - gc1 == 1) { + tranx.data[gc1] = globalTrans->transmissibility(c1, c2); + } + + if (gc2 - gc1 == cartDims[0]) { + trany.data[gc1] = globalTrans->transmissibility(c1, c2); + } + + if (gc2 - gc1 == cartDims[0]*cartDims[1]) { + tranz.data[gc1] = globalTrans->transmissibility(c1, c2); + } + } + } + + return {{"TRANX" , tranx}, + {"TRANY" , trany} , + {"TRANZ" , tranz}}; + } + + Opm::NNC exportNncStructure_() const + { + Opm::NNC nnc = eclState().getInputNNC(); + int nx = eclState().getInputGrid().getNX(); + int ny = eclState().getInputGrid().getNY(); + + const auto& globalGridView = globalGrid_.leafGridView(); + typedef typename Grid::LeafGridView GridView; + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; + ElementMapper globalElemMapper(globalGridView); + + const auto* globalTrans = &(simulator_.gridManager().globalTransmissibility()); + if (!collectToIORank_.isParallel()) { + // in the sequential case we must use the transmissibilites defined by + // the problem. (because in the sequential case, the grid manager does + // not compute "global" transmissibilities for performance reasons. in + // the parallel case, the problem's transmissibilities can't be used + // because this object refers to the distributed grid and we need the + // sequential version here.) + globalTrans = &simulator_.problem().eclTransmissibilities(); + } + + auto elemIt = globalGridView.template begin(); + const auto& elemEndIt = globalGridView.template end(); + for (; elemIt != elemEndIt; ++ elemIt) { + const auto& elem = *elemIt; + + auto isIt = globalGridView.ibegin(elem); + const auto& isEndIt = globalGridView.iend(elem); + for (; isIt != isEndIt; ++ isIt) { + const auto& is = *isIt; + + if (!is.neighbor()) + { + continue; // intersection is on the domain boundary + } + + unsigned c1 = globalElemMapper.index(is.inside()); + unsigned c2 = globalElemMapper.index(is.outside()); + + if (c1 > c2) + { + continue; // we only need to handle each connection once, thank you. + } + + // TODO (?): use the cartesian index mapper to make this code work + // with grids other than Dune::CpGrid. The problem is that we need + // the a mapper for the sequential grid, not for the distributed one. + int cc1 = globalGrid_.globalCell()[c1]; + int cc2 = globalGrid_.globalCell()[c2]; + + if (std::abs(cc1 - cc2) != 1 && + std::abs(cc1 - cc2) != nx && + std::abs(cc1 - cc2) != nx*ny) + { + nnc.addNNC(cc1, cc2, globalTrans->transmissibility(c1, c2)); + } + } + } + return nnc; + } + + const Opm::EclipseState& eclState() const + { return simulator_.gridManager().eclState(); } + const Simulator& simulator_; EclOutputBlackOilModule eclOutputModule_; CollectDataToIORankType collectToIORank_; std::unique_ptr eclIO_; + Grid globalGrid_; }; } // namespace Ewoms