diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index 8082fe2b5..582cd59ef 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -210,7 +210,7 @@ public: simulator_.vanguard().setupTime(); const auto localWellData = simulator_.problem().wellModel().wellData(); - const auto localWBP = data::WellBlockAveragePressures{}; + const auto localWBP = simulator_.problem().wellModel().wellBlockAveragePressures(); const auto localGroupAndNetworkData = simulator_.problem().wellModel() .groupAndNetworkData(reportStepNum); diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 981f66ac7..68d9ea09b 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -38,35 +38,39 @@ #include #include +#include +#include #include #include -#include -#include -#include #include #include + #include #include +#include #include -#include #include #include -#include +#include +#include +#include +#include #include -#include -#include -#include -#include #include #include -#include #include -#include +#include +#include +#include #include +#include #include -#include +#include + +#include #include + #include #include #include @@ -266,6 +270,11 @@ namespace Opm { return wsrpt; } + data::WellBlockAveragePressures wellBlockAveragePressures() const + { + return this->computeWellBlockAveragePressures(); + } + // subtract Binv(D)rw from r; void apply( BVector& r) const; @@ -387,6 +396,13 @@ namespace Opm { std::unique_ptr rateConverter_{}; std::map> regionalAveragePressureCalculator_{}; + struct WBPCalcID + { + std::optional::size_type> openWellIdx_{}; + std::size_t wbpCalcIdx_{}; + }; + + std::vector wbpCalcMap_{}; SimulatorReportSingle last_report_{}; @@ -447,6 +463,16 @@ namespace Opm { // setting the well_solutions_ based on well_state. void updatePrimaryVariables(DeferredLogger& deferred_logger); + void initializeWBPCalculationService(); + + data::WellBlockAveragePressures + computeWellBlockAveragePressures() const; + + ParallelWBPCalculation::EvaluatorFactory + makeWellSourceEvaluatorFactory(const std::vector::size_type wellIdx) const; + + void registerOpenWellsForWBPCalculation(); + void updateAverageFormationFactor(); void computePotentials(const std::size_t widx, diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index fa7fbcc2b..2ba88720e 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -49,6 +49,7 @@ #include #include #include +#include #include #include #include @@ -81,6 +82,7 @@ BlackoilWellModelGeneric(Schedule& schedule, , eclState_(eclState) , comm_(comm) , phase_usage_(phase_usage) + , wbpCalculationService_ { eclState.gridDims(), comm_ } , guideRate_(schedule) , active_wgstate_(phase_usage) , last_valid_wgstate_(phase_usage) @@ -287,57 +289,89 @@ BlackoilWellModelGeneric:: initializeWellPerfData() { well_perf_data_.resize(wells_ecl_.size()); + + this->conn_idx_map_.clear(); + this->conn_idx_map_.reserve(wells_ecl_.size()); + int well_index = 0; for (const auto& well : wells_ecl_) { int connection_index = 0; + // INVALID_ECL_INDEX marks no above perf available int connection_index_above = ParallelWellInfo::INVALID_ECL_INDEX; + well_perf_data_[well_index].clear(); well_perf_data_[well_index].reserve(well.getConnections().size()); - CheckDistributedWellConnections checker(well, local_parallel_well_info_[well_index].get()); + + auto& connIdxMap = this->conn_idx_map_ + .emplace_back(well.getConnections().size()); + + CheckDistributedWellConnections checker { + well, this->local_parallel_well_info_[well_index].get() + }; + bool hasFirstConnection = false; bool firstOpenConnection = true; + auto& parallelWellInfo = this->local_parallel_well_info_[well_index].get(); parallelWellInfo.beginReset(); for (const auto& connection : well.getConnections()) { - const int active_index = compressedIndexForInterior(connection.global_index()); - if (connection.state() == Connection::State::OPEN) { + const auto active_index = + this->compressedIndexForInterior(connection.global_index()); + + const auto connIsOpen = + connection.state() == Connection::State::OPEN; + + if (active_index >= 0) { + connIdxMap.addActiveConnection(connection_index, connIsOpen); + } + + if ((connIsOpen && (active_index >= 0)) || !connIsOpen) { + checker.connectionFound(connection_index); + } + + if (connIsOpen) { if (active_index >= 0) { - if (firstOpenConnection) - { + if (firstOpenConnection) { hasFirstConnection = true; } - checker.connectionFound(connection_index); - PerforationData pd; + + auto pd = PerforationData{}; pd.cell_index = active_index; pd.connection_transmissibility_factor = connection.CF(); pd.satnum_id = connection.satTableId(); pd.ecl_index = connection_index; + well_perf_data_[well_index].push_back(pd); + parallelWellInfo.pushBackEclIndex(connection_index_above, connection_index); } + firstOpenConnection = false; - // Next time this index is the one above as each open connection is - // is stored somehwere. + + // Next time this index is the one above as each open + // connection is stored somewhere. connection_index_above = connection_index; - } else { - checker.connectionFound(connection_index); - if (connection.state() != Connection::State::SHUT) { - OPM_THROW(std::runtime_error, - "Connection state: " + - Connection::State2String(connection.state()) + - " not handled"); - } } - // Note: we rely on the connections being filtered! I.e. there are only connections - // to active cells in the global grid. + else if (connection.state() != Connection::State::SHUT) { + OPM_THROW(std::runtime_error, + fmt::format("Connection state '{}' not handled", + Connection::State2String(connection.state()))); + } + + // Note: we rely on the connections being filtered! I.e., there + // are only connections to active cells in the global grid. ++connection_index; } + parallelWellInfo.endReset(); + checker.checkAllConnectionsFound(); + parallelWellInfo.communicateFirstPerforation(hasFirstConnection); + ++well_index; } } diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp index fa684a43a..dbdc18917 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp @@ -25,20 +25,28 @@ #include -#include #include +#include +#include +#include +#include #include +#include +#include #include #include #include #include +#include #include #include #include +#include #include +#include #include #include #include @@ -419,6 +427,94 @@ protected: std::vector wells_ecl_; std::vector> well_perf_data_; + + /// Connection index mappings + class ConnectionIndexMap + { + public: + /// Constructor. + /// + /// \param[in] numConns Total number of well connections, both open + /// and closed/shut. Typically \code WellConnections::size() \endcode. + explicit ConnectionIndexMap(const std::size_t numConns) + : local_(numConns, -1) + { + this->global_.reserve(numConns); + this->open_.reserve(numConns); + } + + /// Enumerate/map new active connection. + /// + /// \param[in] connIdx Global well connection index. Must be an + /// integer in the range 0..numConns-1. + /// + /// \param[in] connIsOpen Whether or not the connection is + /// open/flowing. + void addActiveConnection(const int connIdx, + const bool connIsOpen) + { + this->local_[connIdx] = + static_cast(this->global_.size()); + + this->global_.push_back(connIdx); + + const auto open_conn_idx = connIsOpen + ? this->num_open_conns_++ + : -1; + + this->open_.push_back(open_conn_idx); + } + + /// Get local connection IDs/indices of every existing well + /// connection. + /// + /// Negative value (-1) for connections that don't intersect the + /// current rank. + const std::vector& local() const + { + return this->local_; + } + + /// Get global connection ID of local (on-rank) connection. + /// + /// \param[in] connIdx Local connection index. + /// + /// \return Global connection ID of \p connIdx. + int global(const int connIdx) const + { + return this->global_[connIdx]; + } + + /// Get open connection ID of local (on-rank) connection. + /// + /// \param[in] connIdx Local connection index. + /// + /// \return Open connection ID of \p connIdx. Integer in the range + /// 0..#open connections - 1 if the connection is open or negative + /// value (-1) otherwise. + int open(const int connIdx) const + { + return this->open_[connIdx]; + } + + private: + /// Local connection IDs/indices of every existing well connection. + /// Negative value (-1) for connections that don't intersect the + /// current rank. + std::vector local_{}; + + /// Global connection index of each on-rank reservoir connection. + /// Reverse/transpose mapping of \c local_. + std::vector global_{}; + + /// Open connection index of each on-rank reservoir connection. + std::vector open_{}; + + /// Number of open connections on this rank. + int num_open_conns_{0}; + }; + + std::vector conn_idx_map_{}; std::function not_on_process_{}; // a vector of all the wells. @@ -430,6 +526,7 @@ protected: std::vector> local_parallel_well_info_; std::vector prod_index_calc_; + mutable ParallelWBPCalculation wbpCalculationService_; std::vector pvt_region_idx_; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 752729edd..abc13524e 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -27,8 +27,11 @@ #include #include #include +#include #include +#include +#include #include #include #include @@ -73,6 +76,48 @@ namespace Opm { this->alternative_well_rate_init_ = EWOMS_GET_PARAM(TypeTag, bool, AlternativeWellRateInit); + + this->wbpCalculationService_ + .localCellIndex([this](const std::size_t globalIndex) + { return this->compressedIndexForInterior(globalIndex); }) + .evalCellSource([this](const int localCell, + PAvgDynamicSourceData::SourceDataSpan sourceTerms) + { + using Item = PAvgDynamicSourceData::SourceDataSpan::Item; + + const auto* intQuants = this->ebosSimulator_.model() + .cachedIntensiveQuantities(localCell, /*timeIndex = */0); + const auto& fs = intQuants->fluidState(); + + sourceTerms.set(Item::PoreVol, intQuants->porosity().value() * + this->ebosSimulator_.model().dofTotalVolume(localCell)); + + constexpr auto io = FluidSystem::oilPhaseIdx; + constexpr auto ig = FluidSystem::gasPhaseIdx; + constexpr auto iw = FluidSystem::waterPhaseIdx; + + // Ideally, these would be 'constexpr'. + const auto haveOil = FluidSystem::phaseIsActive(io); + const auto haveGas = FluidSystem::phaseIsActive(ig); + const auto haveWat = FluidSystem::phaseIsActive(iw); + + auto weightedPhaseDensity = [&fs](const auto ip) + { + return fs.saturation(ip).value() * fs.density(ip).value(); + }; + + if (haveOil) { sourceTerms.set(Item::Pressure, fs.pressure(io).value()); } + else if (haveGas) { sourceTerms.set(Item::Pressure, fs.pressure(ig).value()); } + else { sourceTerms.set(Item::Pressure, fs.pressure(iw).value()); } + + // Strictly speaking, assumes SUM(s[p]) == 1. + auto rho = 0.0; + if (haveOil) { rho += weightedPhaseDensity(io); } + if (haveGas) { rho += weightedPhaseDensity(ig); } + if (haveWat) { rho += weightedPhaseDensity(iw); } + + sourceTerms.set(Item::MixtureDensity, rho); + }); } template @@ -209,6 +254,7 @@ namespace Opm { // We must therefore provide it with updated cell pressures this->initializeWellPerfData(); this->initializeWellState(timeStepIdx, summaryState); + this->initializeWBPCalculationService(); // handling MS well related if (param_.use_multisegment_well_&& anyMSWellOpenLocal()) { // if we use MultisegmentWell model @@ -230,11 +276,16 @@ namespace Opm { } { - const auto& sched_state = this->schedule()[timeStepIdx]; // update VFP properties + const auto& sched_state = this->schedule()[timeStepIdx]; vfp_properties_ = std::make_unique(sched_state.vfpinj(), sched_state.vfpprod(), this->wellState()); + } + + { + const auto& sched_state = this->schedule()[timeStepIdx]; + this->initializeWellProdIndCalculators(); if (sched_state.events().hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX)) { this->runWellPIScaling(timeStepIdx, local_deferredLogger); @@ -799,6 +850,8 @@ namespace Opm { } } } + + this->registerOpenWellsForWBPCalculation(); } @@ -1801,6 +1854,154 @@ namespace Opm { + + + template + void + BlackoilWellModel:: + initializeWBPCalculationService() + { + this->wbpCalcMap_.clear(); + this->wbpCalcMap_.resize(this->wells_ecl_.size()); + + this->registerOpenWellsForWBPCalculation(); + + auto wellID = std::size_t{0}; + for (const auto& well : this->wells_ecl_) { + this->wbpCalcMap_[wellID].wbpCalcIdx_ = this->wbpCalculationService_ + .createCalculator(well, + this->local_parallel_well_info_[wellID], + this->conn_idx_map_[wellID].local(), + this->makeWellSourceEvaluatorFactory(wellID)); + + ++wellID; + } + + this->wbpCalculationService_.defineCommunication(); + } + + + + + + template + data::WellBlockAveragePressures + BlackoilWellModel:: + computeWellBlockAveragePressures() const + { + auto wbpResult = data::WellBlockAveragePressures{}; + + using Calculated = PAvgCalculator::Result::WBPMode; + using Output = data::WellBlockAvgPress::Quantity; + + this->wbpCalculationService_.collectDynamicValues(); + + const auto numWells = this->wells_ecl_.size(); + for (auto wellID = 0*numWells; wellID < numWells; ++wellID) { + const auto calcIdx = this->wbpCalcMap_[wellID].wbpCalcIdx_; + const auto& well = this->wells_ecl_[wellID]; + + if (! well.hasRefDepth()) { + // Can't perform depth correction without at least a + // fall-back datum depth. + continue; + } + + this->wbpCalculationService_ + .inferBlockAveragePressures(calcIdx, well.pavg(), + this->gravity_, + well.getWPaveRefDepth()); + + const auto& result = this->wbpCalculationService_ + .averagePressures(calcIdx); + + auto& reported = wbpResult.values[well.name()]; + + reported[Output::WBP] = result.value(Calculated::WBP); + reported[Output::WBP4] = result.value(Calculated::WBP4); + reported[Output::WBP5] = result.value(Calculated::WBP5); + reported[Output::WBP9] = result.value(Calculated::WBP9); + } + + return wbpResult; + } + + + + + + template + ParallelWBPCalculation::EvaluatorFactory + BlackoilWellModel:: + makeWellSourceEvaluatorFactory(const std::vector::size_type wellIdx) const + { + using Span = PAvgDynamicSourceData::SourceDataSpan; + using Item = typename Span::Item; + + return [wellIdx, this]() -> ParallelWBPCalculation::Evaluator + { + if (! this->wbpCalcMap_[wellIdx].openWellIdx_.has_value()) { + // Well is stopped/shut. Return evaluator for stopped wells. + return []([[maybe_unused]] const int connIdx, Span sourceTerm) + { + // Well/connection is stopped/shut. Set all items to + // zero. + + sourceTerm + .set(Item::Pressure , 0.0) + .set(Item::PoreVol , 0.0) + .set(Item::MixtureDensity, 0.0); + }; + } + + // Well is open. Return an evaluator for open wells/open connections. + return [this, wellPtr = this->well_container_[*this->wbpCalcMap_[wellIdx].openWellIdx_].get()] + (const int connIdx, Span sourceTerm) + { + // Note: The only item which actually matters for the WBP + // calculation at the well reservoir connection level is the + // mixture density. Set other items to zero. + + const auto& connIdxMap = + this->conn_idx_map_[wellPtr->indexOfWell()]; + + const auto rho = wellPtr-> + connectionDensity(connIdxMap.global(connIdx), + connIdxMap.open(connIdx)); + + sourceTerm + .set(Item::Pressure , 0.0) + .set(Item::PoreVol , 0.0) + .set(Item::MixtureDensity, rho); + }; + }; + } + + + + + + template + void + BlackoilWellModel:: + registerOpenWellsForWBPCalculation() + { + assert (this->wbpCalcMap_.size() == this->wells_ecl_.size()); + + for (auto& wbpCalc : this->wbpCalcMap_) { + wbpCalc.openWellIdx_.reset(); + } + + auto openWellIdx = typename std::vector::size_type{0}; + for (const auto* openWell : this->well_container_generic_) { + this->wbpCalcMap_[openWell->indexOfWell()].openWellIdx_ = openWellIdx++; + } + } + + + + + template void BlackoilWellModel:: diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 3a4646462..d3aacc542 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -137,6 +137,9 @@ namespace Opm WellState& well_state, DeferredLogger& deferred_logger) const override; + double connectionDensity(const int globalConnIdx, + const int openConnIdx) const override; + void addWellContributions(SparseMatrixAdapter& jacobian) const override; void addWellPressureEquations(PressureMatrix& mat, diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 73af55623..6a650769d 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -21,8 +21,12 @@ #include #include +#include #include +#include +#include #include + #include #include @@ -703,6 +707,28 @@ namespace Opm + template + double + MultisegmentWell:: + connectionDensity(const int globalConnIdx, + [[maybe_unused]] const int openConnIdx) const + { + // Simple approximation: Mixture density at reservoir connection is + // mixture density at connection's segment. + + const auto segNum = this->wellEcl() + .getConnections()[globalConnIdx].segment(); + + const auto segIdx = this->wellEcl() + .getSegments().segmentNumberToIndex(segNum); + + return this->segments_.density(segIdx).value(); + } + + + + + template void MultisegmentWell:: diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index d80519ef0..6eddd139d 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -185,6 +185,9 @@ namespace Opm WellState& well_state, DeferredLogger& deferred_logger) const override; + virtual double connectionDensity(const int globalConnIdx, + const int openConnIdx) const override; + virtual void addWellContributions(SparseMatrixAdapter& mat) const override; virtual void addWellPressureEquations(PressureMatrix& mat, diff --git a/opm/simulators/wells/StandardWellConnections.hpp b/opm/simulators/wells/StandardWellConnections.hpp index 073e85d13..7c3755d06 100644 --- a/opm/simulators/wells/StandardWellConnections.hpp +++ b/opm/simulators/wells/StandardWellConnections.hpp @@ -73,7 +73,19 @@ public: //! \brief Returns density for first perforation. Scalar rho() const { - return this->perf_densities_.empty() ? 0.0 : perf_densities_[0]; + return this->rho(0); + } + + //! \brief Returns density for specific perforation/connection. + //! + //! \param[in] i Connection index + //! + //! \return Mixture density at connection \p i. + Scalar rho(const typename std::vector::size_type i) const + { + return (i < this->perf_densities_.size()) + ? this->perf_densities_[i] + : 0.0; } //! \brief Returns pressure drop for a given perforation. diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 07ed13a4c..281392ccf 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -1570,6 +1570,23 @@ namespace Opm + + + template + double + StandardWell:: + connectionDensity([[maybe_unused]] const int globalConnIdx, + const int openConnIdx) const + { + return (openConnIdx < 0) + ? 0.0 + : this->connections_.rho(openConnIdx); + } + + + + + template void StandardWell:: diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index ef76195e8..61ad20b62 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -254,6 +254,9 @@ public: WellState& well_state, DeferredLogger& deferred_logger) const = 0; + virtual double connectionDensity(const int globalConnIdx, + const int openConnIdx) const = 0; + /// \brief Wether the Jacobian will also have well contributions in it. virtual bool jacobianContainsWellContributions() const {