diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index f92b2d7cd..7da8f2dec 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -53,6 +53,7 @@ #include #include #include +#include #include #include #include @@ -268,6 +269,7 @@ namespace Opm { std::vector< Well > wells_ecl_; std::vector< std::vector > well_perf_data_; + std::vector< WellProdIndexCalculator > prod_index_calc_; bool wells_active_; @@ -281,6 +283,7 @@ namespace Opm { std::function is_shut_or_defunct_; + void initializeWellProdIndCalculators(); void initializeWellPerfData(); // create the well container @@ -377,6 +380,8 @@ namespace Opm { void calculateEfficiencyFactors(const int reportStepIdx); + void calculateProductivityIndexValues(DeferredLogger& deferred_logger); + // it should be able to go to prepareTimeStep(), however, the updateWellControls() and initPrimaryVariablesEvaluation() // makes it a little more difficult. unless we introduce if (iterationIdx != 0) to avoid doing the above functions // twice at the beginning of the time step diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 859988559..2bc298e80 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -237,6 +237,8 @@ namespace Opm { int globalNumWells = 0; // Make wells_ecl_ contain only this partition's non-shut wells. wells_ecl_ = getLocalNonshutWells(timeStepIdx, globalNumWells); + + this->initializeWellProdIndCalculators(); initializeWellPerfData(); // Wells are active if they are active wells on at least @@ -507,6 +509,8 @@ namespace Opm { const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); checkGconsaleLimits(fieldGroup, well_state_, local_deferredLogger); + this->calculateProductivityIndexValues(local_deferredLogger); + previous_well_state_ = well_state_; Opm::DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger); @@ -566,6 +570,7 @@ namespace Opm { // Make wells_ecl_ contain only this partition's non-shut wells. wells_ecl_ = getLocalNonshutWells(report_step, globalNumWells); + this->initializeWellProdIndCalculators(); initializeWellPerfData(); const int nw = wells_ecl_.size(); @@ -586,6 +591,22 @@ namespace Opm { + template + void + BlackoilWellModel:: + initializeWellProdIndCalculators() + { + this->prod_index_calc_.clear(); + this->prod_index_calc_.reserve(this->wells_ecl_.size()); + for (const auto& well : this->wells_ecl_) { + this->prod_index_calc_.emplace_back(well); + } + } + + + + + template void BlackoilWellModel:: @@ -1439,6 +1460,31 @@ namespace Opm { + + + template + void + BlackoilWellModel:: + calculateProductivityIndexValues(DeferredLogger& deferred_logger) + { + if (! this->localWellsActive()) { + return; + } + + auto piCalc = this->prod_index_calc_.begin(); + for (const auto& wellPtr : this->well_container_) { + wellPtr->updateProductivityIndex(this->ebosSimulator_, + *piCalc, + this->well_state_, + deferred_logger); + ++piCalc; + } + } + + + + + template void BlackoilWellModel:: diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 0cf83fac8..76f44e016 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -44,7 +44,6 @@ namespace Opm using typename Base::RateConverterType; using typename Base::SparseMatrixAdapter; - /// the number of reservior equations using Base::numEq; @@ -168,6 +167,11 @@ namespace Opm const WellState& well_state, Opm::DeferredLogger& deferred_logger) override; // should be const? + virtual void updateProductivityIndex(const Simulator& ebosSimulator, + const WellProdIndexCalculator& wellPICalc, + WellState& well_state, + DeferredLogger& deferred_logger) const override; + virtual void addWellContributions(SparseMatrixAdapter& jacobian) const override; /// number of segments for this well diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index b3bb660ce..f45eec18a 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1176,6 +1176,94 @@ namespace Opm + template + void + MultisegmentWell:: + updateProductivityIndex(const Simulator& ebosSimulator, + const WellProdIndexCalculator& wellPICalc, + WellState& well_state, + DeferredLogger&) const + { + auto recipFVF = [&ebosSimulator, this](const int perf, const int p) -> double + { + const auto cell_idx = this->well_cells_[perf]; + const auto& fs = ebosSimulator.model() + .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState(); + + return fs.invB(p).value(); + }; + + auto Rs = [&ebosSimulator, this](const int perf) -> double + { + const auto cell_idx = this->well_cells_[perf]; + const auto& fs = ebosSimulator.model() + .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState(); + + return fs.Rs().value(); + }; + + auto Rv = [&ebosSimulator, this](const int perf) -> double + { + const auto cell_idx = this->well_cells_[perf]; + const auto& fs = ebosSimulator.model() + .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState(); + + return fs.Rv().value(); + }; + + const auto& pu = this->phaseUsage(); + const int np = this->number_of_phases_; + + auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0]; + auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0]; + + for (int p = 0; p < np; ++p) { + wellPI[p] = 0.0; + } + + const auto& allConn = this->well_ecl_.getConnections(); + const auto nPerf = allConn.size(); + for (auto allPerfID = 0*nPerf, subsetPerfID = 0*nPerf; allPerfID < nPerf; ++allPerfID) { + if (allConn[allPerfID].state() == Connection::State::SHUT) { + continue; + } + + std::vector mob(num_components_, 0.0); + getMobility(ebosSimulator, static_cast(subsetPerfID), mob); + + for (int p = 0; p < np; ++p) { + // Note: E100's notion of PI value phase mobility includes + // the reciprocal FVF. + const auto connMob = mob[ flowPhaseToEbosCompIdx(p) ].value() * recipFVF(subsetPerfID, p); + connPI[p] = wellPICalc.connectionProdIndStandard(allPerfID, connMob); + } + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) + { + const auto io = pu.phase_pos[Oil]; + const auto ig = pu.phase_pos[Gas]; + + const auto vapoil = connPI[ig] * Rv(subsetPerfID); + const auto disgas = connPI[io] * Rs(subsetPerfID); + + connPI[io] += vapoil; + connPI[ig] += disgas; + } + + for (int p = 0; p < np; ++p) { + wellPI[p] += connPI[p]; + } + + ++subsetPerfID; + connPI += np; + } + } + + + + + template void MultisegmentWell:: diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index c3aa8d6df..341d3dd95 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -27,9 +27,10 @@ #include #endif +#include #include #include -#include +#include #include #include @@ -224,6 +225,11 @@ namespace Opm const WellState& well_state, Opm::DeferredLogger& deferred_logger) override; // should be const? + virtual void updateProductivityIndex(const Simulator& ebosSimulator, + const WellProdIndexCalculator& wellPICalc, + WellState& well_state, + DeferredLogger& deferred_logger) const override; + virtual void addWellContributions(SparseMatrixAdapter& mat) const override; // iterate well equations with the specified control until converged @@ -315,7 +321,6 @@ namespace Opm using Base::wpolymer; using Base::wfoam; using Base::scalingFactor; - using Base::scaleProductivityIndex; using Base::mostStrictBhpFromBhpLimits; // protected member variables from the Base class diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 43752fc18..01986fa94 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -602,9 +602,6 @@ namespace Opm well_state.wellDissolvedGasRates()[index_of_well_] = 0.; const int np = number_of_phases_; - for (int p = 0; p < np; ++p) { - well_state.productivityIndex()[np*index_of_well_ + p] = 0.; - } std::vector connectionRates = connectionRates_; // Copy to get right size. for (int perf = 0; perf < number_of_perforations_; ++perf) { @@ -683,8 +680,6 @@ namespace Opm } catch( ... ) { OPM_DEFLOG_THROW(Opm::NumericalIssue,"Error when inverting local well equations for well " + name(), deferred_logger); } - - } @@ -703,7 +698,6 @@ namespace Opm Opm::DeferredLogger& deferred_logger) const { const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator); - const int np = number_of_phases_; const EvalWell& bhp = getBhp(); const int cell_idx = well_cells_[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); @@ -860,27 +854,6 @@ namespace Opm // Store the perforation pressure for later usage. well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf]; - - // Compute Productivity index if asked for - const auto& pu = phaseUsage(); - const Opm::SummaryConfig& summaryConfig = ebosSimulator.vanguard().summaryConfig(); - const Opm::Schedule& schedule = ebosSimulator.vanguard().schedule(); - for (int p = 0; p < np; ++p) { - if ( (pu.phase_pos[Water] == p && (summaryConfig.hasSummaryKey("WPIW:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name()))) - || (pu.phase_pos[Oil] == p && (summaryConfig.hasSummaryKey("WPIO:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name()))) - || (pu.phase_pos[Gas] == p && summaryConfig.hasSummaryKey("WPIG:" + name()))) { - - const unsigned int compIdx = flowPhaseToEbosCompIdx(p); - const auto& fs = intQuants.fluidState(); - Eval perf_pressure = getPerfCellPressure(fs); - const double drawdown = well_state.perfPress()[first_perf_ + perf] - perf_pressure.value(); - const bool new_well = schedule.hasWellGroupEvent(name(), ScheduleEvents::NEW_WELL, current_step_); - double productivity_index = cq_s[compIdx].value() / drawdown; - scaleProductivityIndex(perf, productivity_index, new_well, deferred_logger); - well_state.productivityIndex()[np*index_of_well_ + p] += productivity_index; - } - } - } @@ -2307,6 +2280,94 @@ namespace Opm + template + void + StandardWell:: + updateProductivityIndex(const Simulator& ebosSimulator, + const WellProdIndexCalculator& wellPICalc, + WellState& well_state, + DeferredLogger& deferred_logger) const + { + auto recipFVF = [&ebosSimulator, this](const int perf, const int p) -> double + { + const auto cell_idx = this->well_cells_[perf]; + const auto& fs = ebosSimulator.model() + .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState(); + + return fs.invB(p).value(); + }; + + auto Rs = [&ebosSimulator, this](const int perf) -> double + { + const auto cell_idx = this->well_cells_[perf]; + const auto& fs = ebosSimulator.model() + .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState(); + + return fs.Rs().value(); + }; + + auto Rv = [&ebosSimulator, this](const int perf) -> double + { + const auto cell_idx = this->well_cells_[perf]; + const auto& fs = ebosSimulator.model() + .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState(); + + return fs.Rv().value(); + }; + + const auto& pu = this->phaseUsage(); + const int np = this->number_of_phases_; + + auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0]; + auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0]; + + for (int p = 0; p < np; ++p) { + wellPI[p] = 0.0; + } + + const auto& allConn = this->well_ecl_.getConnections(); + const auto nPerf = allConn.size(); + for (auto allPerfID = 0*nPerf, subsetPerfID = 0*nPerf; allPerfID < nPerf; ++allPerfID) { + if (allConn[allPerfID].state() == Connection::State::SHUT) { + continue; + } + + std::vector mob(num_components_, {numWellEq_ + numEq, 0.0}); + getMobility(ebosSimulator, static_cast(subsetPerfID), mob, deferred_logger); + + for (int p = 0; p < np; ++p) { + // Note: E100's notion of PI value phase mobility includes + // the reciprocal FVF. + const auto connMob = mob[ flowPhaseToEbosCompIdx(p) ].value() * recipFVF(subsetPerfID, p); + connPI[p] = wellPICalc.connectionProdIndStandard(allPerfID, connMob); + } + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) + { + const auto io = pu.phase_pos[Oil]; + const auto ig = pu.phase_pos[Gas]; + + const auto vapoil = connPI[ig] * Rv(subsetPerfID); + const auto disgas = connPI[io] * Rs(subsetPerfID); + + connPI[io] += vapoil; + connPI[ig] += disgas; + } + + for (int p = 0; p < np; ++p) { + wellPI[p] += connPI[p]; + } + + ++subsetPerfID; + connPI += np; + } + } + + + + + template void StandardWell:: diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 4ca828540..f90ca706d 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -39,6 +39,7 @@ #include #include #include +#include #include #include @@ -221,6 +222,11 @@ namespace Opm const WellState& well_state, Opm::DeferredLogger& deferred_logger) = 0; // should be const? + virtual void updateProductivityIndex(const Simulator& ebosSimulator, + const WellProdIndexCalculator& wellPICalc, + WellState& well_state, + DeferredLogger& deferred_logger) const = 0; + /// \brief Wether the Jacobian will also have well contributions in it. virtual bool jacobianContainsWellContributions() const { @@ -513,14 +519,8 @@ namespace Opm const std::vector& B_avg, Opm::DeferredLogger& deferred_logger); - void scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger) const; - void initCompletions(); - // count the number of times an output log message is created in the productivity - // index calculations - mutable int well_productivity_index_logger_counter_; - bool checkConstraints(WellState& well_state, const Schedule& schedule, const SummaryState& summaryState, diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 675e5afbf..91ee4dd91 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -85,8 +85,6 @@ namespace Opm connectionRates_.resize(number_of_perforations_); - well_productivity_index_logger_counter_ = 0; - wellIsStopped_ = false; if (well.getStatus() == Well::Status::STOP) { wellIsStopped_ = true; @@ -1380,40 +1378,6 @@ namespace Opm } } - template - void - WellInterface::scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger) const - { - const auto& connection = well_ecl_.getConnections()[(*perf_data_)[perfIdx].ecl_index]; - if (well_ecl_.getDrainageRadius() < 0) { - if (new_well && perfIdx == 0) { - deferred_logger.warning("PRODUCTIVITY_INDEX_WARNING", "Negative drainage radius not supported. The productivity index is set to zero"); - } - productivity_index = 0.0; - return; - } - - if (connection.r0() > well_ecl_.getDrainageRadius()) { - if (new_well && well_productivity_index_logger_counter_ < 1) { - deferred_logger.info("PRODUCTIVITY_INDEX_INFO", "The effective radius is larger than the well drainage radius for well " + name() + - " They are set to equal in the well productivity index calculations"); - well_productivity_index_logger_counter_++; - } - return; - } - - // For zero drainage radius the productivity index is just the transmissibility times the mobility - if (well_ecl_.getDrainageRadius() == 0) { - return; - } - - // Scale the productivity index to account for the drainage radius. - // Assumes steady radial flow only valied for horizontal wells - productivity_index *= - (std::log(connection.r0() / connection.rw()) + connection.skinFactor()) / - (std::log(well_ecl_.getDrainageRadius() / connection.rw()) + connection.skinFactor()); - } - template void WellInterface::addCellRates(RateVector& rates, int cellIdx) const diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp index 1e73cbced..10c0ec25f 100644 --- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp +++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp @@ -29,15 +29,15 @@ #include -#include -#include -#include -#include -#include #include #include +#include #include +#include #include +#include +#include +#include namespace Opm { @@ -161,6 +161,7 @@ namespace Opm perfRateSolvent_.clear(); perfRateSolvent_.resize(nperf, 0.0); productivity_index_.resize(nw * np, 0.0); + conn_productivity_index_.resize(nperf * np, 0.0); well_potentials_.resize(nw * np, 0.0); perfRatePolymer_.clear(); @@ -515,16 +516,20 @@ namespace Opm using rt = data::Rates::opt; std::vector< rt > phs( np ); + std::vector pi(np); if( pu.phase_used[Water] ) { phs.at( pu.phase_pos[Water] ) = rt::wat; + pi .at( pu.phase_pos[Water] ) = rt::productivity_index_water; } if( pu.phase_used[Oil] ) { phs.at( pu.phase_pos[Oil] ) = rt::oil; + pi .at( pu.phase_pos[Oil] ) = rt::productivity_index_oil; } if( pu.phase_used[Gas] ) { phs.at( pu.phase_pos[Gas] ) = rt::gas; + pi .at( pu.phase_pos[Gas] ) = rt::productivity_index_gas; } /* this is a reference or example on **how** to convert from @@ -612,13 +617,14 @@ namespace Opm size_t local_comp_index = 0; for( auto& comp : well.connections) { - const auto rates = this->perfPhaseRates().begin() - + (np * wt.second[ 1 ]) - + (np * local_comp_index); + const auto connPhaseOffset = np * (wt.second[1] + local_comp_index); + const auto rates = this->perfPhaseRates().begin() + connPhaseOffset; + const auto connPI = this->connectionProductivityIndex().begin() + connPhaseOffset; for( int i = 0; i < np; ++i ) { - comp.rates.set( phs[ i ], *(rates + i) ); + comp.rates.set( phs[ i ], *(rates + i) ); + comp.rates.set( pi [ i ], *(connPI + i) ); } if ( pu.has_polymer ) { comp.rates.set( rt::polymer, this->perfRatePolymer()[wt.second[1] + local_comp_index]); @@ -981,6 +987,14 @@ namespace Opm return productivity_index_; } + std::vector& connectionProductivityIndex() { + return this->conn_productivity_index_; + } + + const std::vector& connectionProductivityIndex() const { + return this->conn_productivity_index_; + } + std::vector& wellPotentials() { return well_potentials_; } @@ -1262,6 +1276,9 @@ namespace Opm // Productivity Index std::vector productivity_index_; + // Connection-level Productivity Index + std::vector conn_productivity_index_; + // Well potentials std::vector well_potentials_;