diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp index 9c877c791..7d3d704c2 100644 --- a/opm/autodiff/BlackoilWellModel_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -396,6 +396,10 @@ namespace Opm { initial_step_ = false; } + + + + template std::vector::WellInterfacePtr > BlackoilWellModel:: @@ -429,21 +433,33 @@ namespace Opm { const Well* well_ecl = wells_ecl_[index_well]; - // well is closed due to economical reasons - if (wellTestState_.hasWell(well_name, WellTestConfig::Reason::ECONOMIC)) { - if( well_ecl->getAutomaticShutIn() ) { - // shut wells are not added to the well container - well_state_.bhp()[w] = 0; + // A new WCON keywords can re-open a well that was closed/shut due to Physical limit + if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL ) ) { + // TODO: more checking here, to make sure this standard more specific and complete + // maybe there is some WCON keywords will not open the well + if (well_state_.effectiveEventsOccurred(w) ) { + wellTestState_.openWell(well_name); + } + } + + // TODO: should we do this for all kinds of closing reasons? + // something like wellTestState_.hasWell(well_name)? + if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::ECONOMIC) || + wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL) ) { + if( well_ecl->getAutomaticShutIn() ) { + // shut wells are not added to the well container + // TODO: make a function from well_state side to handle the following + well_state_.thp()[w] = 0.; + well_state_.bhp()[w] = 0.; const int np = numPhases(); for (int p = 0; p < np; ++p) { - well_state_.wellRates()[np * w + p] = 0; - } - continue; - } - else { - // close wells are added to the container but marked as closed - struct WellControls* well_controls = wells()->ctrls[w]; - well_controls_stop_well(well_controls); + well_state_.wellRates()[np * w + p] = 0.; + } + continue; + } else { + // close wells are added to the container but marked as closed + struct WellControls* well_controls = wells()->ctrls[w]; + well_controls_stop_well(well_controls); } } @@ -803,7 +819,9 @@ namespace Opm { // Get global (from all processes) convergence report. ConvergenceReport local_report; for (const auto& well : well_container_) { - local_report += well->getWellConvergence(B_avg); + if (well->isOperable() ) { + local_report += well->getWellConvergence(B_avg); + } } ConvergenceReport report = gatherConvergenceReport(local_report); @@ -828,9 +846,10 @@ namespace Opm { BlackoilWellModel:: calculateExplicitQuantities() const { - for (auto& well : well_container_) { - well->calculateExplicitQuantities(ebosSimulator_, well_state_); - } + // TODO: checking isOperable() ? + for (auto& well : well_container_) { + well->calculateExplicitQuantities(ebosSimulator_, well_state_); + } } @@ -934,6 +953,10 @@ namespace Opm { // process group control related prepareGroupControl(); + for (const auto& well : well_container_) { + well->checkWellOperability(ebosSimulator_, well_state_); + } + // since the controls are all updated, we should update well_state accordingly for (const auto& well : well_container_) { const int w = well->indexOfWell(); @@ -941,16 +964,22 @@ namespace Opm { const int control = well_controls_get_current(wc); well_state_.currentControls()[w] = control; + if (!well->isOperable() ) continue; + if (well_state_.effectiveEventsOccurred(w) ) { well->updateWellStateWithTarget(ebosSimulator_, well_state_); } // there is no new well control change input within a report step, // so next time step, the well does not consider to have effective events anymore + // TODO: if we can know whether this is the first time step within the report step, + // we do not need to set it to false + // TODO: we should do this at the end of the time step in case we will need it within + // this time step somewhere if (well_state_.effectiveEventsOccurred(w) ) { well_state_.setEffectiveEventsOccurred(w, false); } - } // end of for (int w = 0; w < nw; ++w) + } // end of for (const auto& well : well_container_) updatePrimaryVariables(); } diff --git a/opm/autodiff/MultisegmentWell.hpp b/opm/autodiff/MultisegmentWell.hpp index 40496a55e..318d06609 100644 --- a/opm/autodiff/MultisegmentWell.hpp +++ b/opm/autodiff/MultisegmentWell.hpp @@ -330,6 +330,11 @@ namespace Opm // handling the overshooting and undershooting of the fractions void processFractions(const int seg) const; + // checking the operability of the well based on current reservoir condition + // it is not implemented for multisegment well yet + virtual void checkWellOperability(const Simulator& ebos_simulator, + const WellState& well_state) override; + void updateWellStateFromPrimaryVariables(WellState& well_state) const; bool frictionalPressureLossConsidered() const; diff --git a/opm/autodiff/MultisegmentWell_impl.hpp b/opm/autodiff/MultisegmentWell_impl.hpp index 98271352c..fd3a1362e 100644 --- a/opm/autodiff/MultisegmentWell_impl.hpp +++ b/opm/autodiff/MultisegmentWell_impl.hpp @@ -1642,6 +1642,21 @@ namespace Opm + template + void + MultisegmentWell:: + checkWellOperability(const Simulator& /* ebos_simulator */, + const WellState& /* well_state */) + { + const std::string msg = "Support of well operatability checking for mutlisegment wells is not implemented " + "yet, checkWellOperability() for " + name() + " will do nothing"; + OpmLog::warning("NO_OPERATABILITY_CHECKING_MS_WELLS", msg); + } + + + + + template void MultisegmentWell:: diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index a2d2b27ff..b5ff47380 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -264,8 +264,6 @@ namespace Opm EvalWell extendEval(const Eval& in) const; - bool crossFlowAllowed(const Simulator& ebosSimulator) const; - // xw = inv(D)*(rw - C*x) void recoverSolutionWell(const BVector& x, BVectorWell& xw) const; @@ -350,6 +348,17 @@ namespace Opm // updating the inflow based on the current reservoir condition void updateIPR(const Simulator& ebos_simulator) const; + // check whether the well is operable under the current reservoir condition + // mostly related to BHP limit and THP limit + virtual void checkWellOperability(const Simulator& ebos_simulator, + const WellState& well_state) override; + + // check whether the well is operable under BHP limit with current reservoir condition + void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator); + + // check whether the well is operable under THP limit with current reservoir condition + void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator); + // update WellState based on IPR and associated VFP table void updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator, WellState& well_state) const; @@ -357,6 +366,21 @@ namespace Opm void updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator, WellState& well_state) const; + // for a well, when all drawdown are in the wrong direction, then this well will not + // be able to produce/inject . + bool allDrawDownWrongDirection(const Simulator& ebos_simulator) const; + + // whether the well can produce / inject based on the current well state (bhp) + bool canProduceInjectWithCurrentBhp(const Simulator& ebos_simulator, + const WellState& well_state); + + // turn on crossflow to avoid singular well equations + // when the well is banned from cross-flow and the BHP is not properly initialized, + // we turn on crossflow to avoid singular well equations. It can result in wrong-signed + // well rates, it can cause problem for THP calculation + // TODO: looking for better alternative to avoid wrong-signed well rates + bool openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const; + // calculate the BHP from THP target based on IPR // TODO: we need to check the operablility here first, if not operable, then maybe there is // no point to do this diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 5da6bf514..5d944b626 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -440,9 +440,9 @@ namespace Opm WellState& well_state) { - const Opm::SummaryConfig& summaryConfig = ebosSimulator.vanguard().summaryConfig(); + checkWellOperability(ebosSimulator, well_state); - const int np = number_of_phases_; + if (!this->isOperable()) return; // clear all entries duneB_ = 0.0; @@ -453,7 +453,7 @@ namespace Opm // TODO: it probably can be static member for StandardWell const double volume = 0.002831684659200; // 0.1 cu ft; - const bool allow_cf = crossFlowAllowed(ebosSimulator); + const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator); const EvalWell& bhp = getBhp(); @@ -461,6 +461,7 @@ namespace Opm well_state.wellVaporizedOilRates()[index_of_well_] = 0.; 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.; } @@ -588,7 +589,8 @@ namespace Opm 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 auto& pu = phaseUsage(); + const Opm::SummaryConfig& summaryConfig = ebosSimulator.vanguard().summaryConfig(); 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()))) @@ -747,45 +749,6 @@ namespace Opm - template - bool - StandardWell:: - crossFlowAllowed(const Simulator& ebosSimulator) const - { - if (getAllowCrossFlow()) { - return true; - } - - // TODO: investigate the justification of the following situation - - // check for special case where all perforations have cross flow - // then the wells must allow for cross flow - for (int perf = 0; perf < number_of_perforations_; ++perf) { - const int cell_idx = well_cells_[perf]; - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - const auto& fs = intQuants.fluidState(); - const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); - const EvalWell& bhp = getBhp(); - - // Pressure drawdown (also used to determine direction of flow) - const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf]; - const EvalWell drawdown = pressure - well_pressure; - - if (drawdown.value() < 0 && well_type_ == INJECTOR) { - return false; - } - - if (drawdown.value() > 0 && well_type_ == PRODUCER) { - return false; - } - } - return true; - } - - - - - template void StandardWell:: @@ -860,6 +823,8 @@ namespace Opm updateWellState(const BVectorWell& dwells, WellState& well_state) const { + if (!this->isOperable()) return; + updatePrimaryVariablesNewton(dwells, well_state); updateWellStateFromPrimaryVariables(well_state); @@ -1144,13 +1109,21 @@ namespace Opm // TODO: we should address this in a function updateWellStateWithBHPTarget. // TODO: however, the reason that this one minght not be that critical with // TODO: the effects remaining to be investigated. - break; + case THP: { - // TODO: adding the checking for the operability - // TODO: should we do updateIPR before this or within the related functions - updateIPR(ebos_simulator); - updateWellStateWithTHPTargetIPR(ebos_simulator, well_state); + assert(this->isOperable() ); + + // when a well can not work under THP target, it switches to BHP control + if (this->operability_status_.isOperableUnderTHPLimit() ) { + updateWellStateWithTHPTargetIPR(ebos_simulator, well_state); + } else { // go to BHP limit + assert(this->operability_status_.isOperableUnderBHPLimit() ); + + OpmLog::info("well " + name() + " can not work with THP target, switching to BHP control"); + + well_state.bhp()[well_index] = mostStrictBhpFromBhpLimits(); + } break; } @@ -1267,7 +1240,8 @@ namespace Opm // it should not be negative anyway. If it is negative, we might need to re-formulate // to taking into consideration the crossflow here. if (pressure_diff <= 0.) { - OpmLog::warning("NON_POSITIVE_DRAWDOWN_IPR", "non-positive drawdown found when updateIPR for well " + name()); + OpmLog::warning("NON_POSITIVE_DRAWDOWN_IPR", + "non-positive drawdown found when updateIPR for well " + name()); } // the well index associated with the connection @@ -1316,6 +1290,221 @@ namespace Opm + template + void + StandardWell:: + checkWellOperability(const Simulator& ebos_simulator, + const WellState& well_state) + { + // focusing on PRODUCER for now + if (well_type_ == INJECTOR) { + return; + } + + if (!this->underPredictionMode() ) { + return; + } + + const bool old_well_operable = this->operability_status_.isOperable(); + + this->operability_status_.reset(); + + updateIPR(ebos_simulator); + + // checking the BHP limit related + checkOperabilityUnderBHPLimitProducer(ebos_simulator); + + // checking whether the well can operate under the THP constraints. + if (this->wellHasTHPConstraints()) { + this->operability_status_.has_thp_constaint = true; + checkOperabilityUnderTHPLimitProducer(ebos_simulator); + this->operability_status_.can_produce_inject_with_current_bhp = + canProduceInjectWithCurrentBhp(ebos_simulator, well_state); + } + + const bool well_operable = this->operability_status_.isOperable(); + + if (!well_operable && old_well_operable) { + OpmLog::info(" well " + name() + " gets SHUT during iteration "); + } else if (well_operable && !old_well_operable) { + OpmLog::info(" well " + name() + " gets REVIVED during iteration "); + } + } + + + + + + template + void + StandardWell:: + checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator) + { + const double bhp_limit = mostStrictBhpFromBhpLimits(); + // Crude but works: default is one atmosphere. + // TODO: a better way to detect whether the BHP is defaulted or not + const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa; + if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints() ) { + // if the BHP limit is not defaulted or the well does not have a THP limit + // we need to check the BHP limit + + for (int p = 0; p < number_of_phases_; ++p) { + const double temp = ipr_a_[p] - ipr_b_[p] * bhp_limit; + if (temp < 0.) { + this->operability_status_.operable_under_only_bhp_limit = false; + break; + } + } + + // checking whether running under BHP limit will violate THP limit + if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints()) { + // option 1: calculate well rates based on the BHP limit. + // option 2: stick with the above IPR curve + // we use IPR here + std::vector well_rates_bhp_limit; + computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit); + + const double thp = calculateThpFromBhp(well_rates_bhp_limit, bhp_limit); + const double thp_limit = this->getTHPConstraint(); + + if (thp < thp_limit) { + this->operability_status_.obey_thp_limit_under_bhp_limit = false; + } + } + } else { + // defaulted BHP and there is a THP constraint + // default BHP limit is about 1 atm. + // when applied the hydrostatic pressure correction dp, + // most likely we get a negative value (bhp + dp)to search in the VFP table, + // which is not desirable. + // we assume we can operate under defaulted BHP limit and will violate the THP limit + // when operating under defaulted BHP limit. + this->operability_status_.operable_under_only_bhp_limit = true; + this->operability_status_.obey_thp_limit_under_bhp_limit = false; + } + } + + + + + + template + void + StandardWell:: + checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator) + { + const double obtain_bhp = calculateBHPWithTHPTargetIPR(); + + if (obtain_bhp > 0.) { + this->operability_status_.can_obtain_bhp_with_thp_limit = true; + + const double bhp_limit = mostStrictBhpFromBhpLimits(); + this->operability_status_.obey_bhp_limit_with_thp_limit = (obtain_bhp >= bhp_limit); + + const double thp_limit = this->getTHPConstraint(); + if (obtain_bhp < thp_limit) { + const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(obtain_bhp, unit::barsa)) + + " bars is SMALLER than thp limit " + + std::to_string(unit::convert::to(thp_limit, unit::barsa)) + + " bars as a producer for well " + name(); + OpmLog::debug(msg); + } + } else { + this->operability_status_.can_obtain_bhp_with_thp_limit = false; + const double thp_limit = this->getTHPConstraint(); + OpmLog::debug(" COULD NOT find bhp value under thp_limit " + + std::to_string(unit::convert::to(thp_limit, unit::barsa)) + + " bars for well " + name() + ", the well might need to be closed "); + this->operability_status_.obey_bhp_limit_with_thp_limit = false; + } + } + + + + + + template + bool + StandardWell:: + allDrawDownWrongDirection(const Simulator& ebos_simulator) const + { + bool all_drawdown_wrong_direction = true; + + for (int perf = 0; perf < number_of_perforations_; ++perf) { + const int cell_idx = well_cells_[perf]; + const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + const double pressure = (fs.pressure(FluidSystem::oilPhaseIdx)).value(); + const double bhp = getBhp().value(); + + // Pressure drawdown (also used to determine direction of flow) + const double well_pressure = bhp + perf_pressure_diffs_[perf]; + const double drawdown = pressure - well_pressure; + + // for now, if there is one perforation can produce/inject in the correct + // direction, we consider this well can still produce/inject. + // TODO: it can be more complicated than this to cause wrong-signed rates + if ( (drawdown < 0. && well_type_ == INJECTOR) || + (drawdown > 0. && well_type_ == PRODUCER) ) { + all_drawdown_wrong_direction = false; + break; + } + } + + return all_drawdown_wrong_direction; + } + + + + + template + bool + StandardWell:: + canProduceInjectWithCurrentBhp(const Simulator& ebos_simulator, + const WellState& well_state) + { + const double bhp = well_state.bhp()[index_of_well_]; + std::vector well_rates; + computeWellRatesWithBhp(ebos_simulator, bhp, well_rates); + + const double sign = (well_type_ == PRODUCER) ? -1. : 1.; + const double threshold = sign * std::numeric_limits::min(); + + bool can_produce_inject = false; + for (const auto value : well_rates) { + if (well_type_ == PRODUCER && value < threshold) { + can_produce_inject = true; + break; + } else if (well_type_ == INJECTOR && value > threshold) { + can_produce_inject = true; + break; + } + } + + if (!can_produce_inject) { + OpmLog::debug(" well " + name() + " CANNOT produce or inejct "); + } + + return can_produce_inject; + } + + + + + + template + bool + StandardWell:: + openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const + { + return !getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator); + } + + + + + template void StandardWell:: @@ -1349,6 +1538,8 @@ namespace Opm const double bhp = calculateBHPWithTHPTargetIPR(); + assert(bhp > 0.0); + well_state.bhp()[index_of_well_] = bhp; // TODO: explicit quantities are always tricky for this type of situation @@ -1401,9 +1592,6 @@ namespace Opm const double obtain_bhp = vfp_properties_->getProd()->calculateBhpWithTHPTarget(ipr_a_, ipr_b_, bhp_limit, thp_table_id, thp_target, alq, dp); - // we should have made sure that this well should be operable under THP limit now - assert(obtain_bhp > 0.); - return obtain_bhp; } @@ -1726,7 +1914,7 @@ namespace Opm switch(well_controls_get_current_type(well_controls_)) { case THP: type = CR::WellFailure::Type::ControlTHP; - control_tolerance = 1.e3; // 0.01 bar + control_tolerance = 1.e4; // 0.1 bar break; case BHP: // pressure type of control type = CR::WellFailure::Type::ControlBHP; @@ -1816,6 +2004,8 @@ namespace Opm StandardWell:: solveEqAndUpdateWellState(WellState& well_state) { + if (!this->isOperable()) return; + // We assemble the well equations, then we check the convergence, // which is why we do not put the assembleWellEq here. BVectorWell dx_well(1); @@ -1861,6 +2051,8 @@ namespace Opm StandardWell:: apply(const BVector& x, BVector& Ax) const { + if (!this->isOperable()) return; + if ( param_.matrix_add_well_contributions_ ) { // Contributions are already in the matrix itself @@ -1889,6 +2081,8 @@ namespace Opm StandardWell:: apply(BVector& r) const { + if (!this->isOperable()) return; + assert( invDrw_.size() == invDuneD_.N() ); // invDrw_ = invDuneD_ * resWell_ @@ -1906,6 +2100,8 @@ namespace Opm StandardWell:: recoverSolutionWell(const BVector& x, BVectorWell& xw) const { + if (!this->isOperable()) return; + BVectorWell resWell = resWell_; // resWell = resWell - B * x duneB_.mmv(x, resWell); @@ -1923,6 +2119,8 @@ namespace Opm recoverWellSolutionAndUpdateWellState(const BVector& x, WellState& well_state) const { + if (!this->isOperable()) return; + BVectorWell xw(1); recoverSolutionWell(x, xw); updateWellState(xw, well_state); @@ -1942,7 +2140,7 @@ namespace Opm const int np = number_of_phases_; well_flux.resize(np, 0.0); - const bool allow_cf = crossFlowAllowed(ebosSimulator); + const bool allow_cf = getAllowCrossFlow(); for (int perf = 0; perf < number_of_perforations_; ++perf) { const int cell_idx = well_cells_[perf]; @@ -2127,6 +2325,8 @@ namespace Opm StandardWell:: updatePrimaryVariables(const WellState& well_state) const { + if (!this->isOperable()) return; + const int well_index = index_of_well_; const int np = number_of_phases_; @@ -2337,7 +2537,8 @@ namespace Opm return; } // compute the well water velocity with out shear effects. - const bool allow_cf = crossFlowAllowed(ebos_simulator); + // TODO: do we need to turn on crossflow here? + const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator); const EvalWell& bhp = getBhp(); std::vector cq_s(num_components_,0.0); double perf_dis_gas_rate = 0.; diff --git a/opm/autodiff/VFPProdProperties.cpp b/opm/autodiff/VFPProdProperties.cpp index f86c0be40..89c85adfb 100644 --- a/opm/autodiff/VFPProdProperties.cpp +++ b/opm/autodiff/VFPProdProperties.cpp @@ -229,10 +229,10 @@ calculateBhpWithTHPTarget(const std::vector& ipr_a, detail::RateBhpPair{flo_bhp_limit, bhp_limit} }; double obtain_bhp = 0.; - const bool obtain_solution_with_thp_limit = detail::findIntersectionForBhp(ratebhp_samples, ratebhp_twopoints_ipr, obtain_bhp); + const bool can_obtain_bhp_with_thp_limit = detail::findIntersectionForBhp(ratebhp_samples, ratebhp_twopoints_ipr, obtain_bhp); - // \Note: assuming not that negative BHP does not make sense - if (obtain_solution_with_thp_limit && obtain_bhp > 0.) { + // \Note: assuming that negative BHP does not make sense + if (can_obtain_bhp_with_thp_limit && obtain_bhp > 0.) { // getting too high bhp that might cause negative rates (rates in the undesired direction) if (obtain_bhp >= bhp_safe_limit) { const std::string msg (" We are getting a too high BHP value from the THP constraint, which may " @@ -252,5 +252,4 @@ calculateBhpWithTHPTarget(const std::vector& ipr_a, } - } diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index c084e9ab0..272809588 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -223,6 +223,10 @@ namespace Opm void updatePerforatedCell(std::vector& is_cell_perforated); + virtual void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state) = 0; + + // whether the well is operable + bool isOperable() const; protected: @@ -350,6 +354,10 @@ namespace Opm // whether a well is specified with a non-zero and valid VFP table number bool isVFPActive() const; + struct OperabilityStatus; + + OperabilityStatus operability_status_; + void wellTestingEconomic(Simulator& simulator, const std::vector& B_avg, const double simulation_time, const int report_step, const bool terminal_output, const WellState& well_state, WellTestState& welltest_state); @@ -359,6 +367,11 @@ namespace Opm const bool write_message_to_opmlog, WellTestState& well_test_state) const; + void updateWellTestStatePhysical(const WellState& well_state, + const double simulation_time, + const bool write_message_to_opmlog, + WellTestState& well_test_state) const; + void solveWellForTesting(Simulator& ebosSimulator, WellState& well_state, const std::vector& B_avg, bool terminal_output); void scaleProductivityIndex(const int perfIdx, double& productivity_index) const; @@ -366,6 +379,62 @@ namespace Opm }; + + + + // definition of the struct OperabilityStatus + template + struct + WellInterface:: + OperabilityStatus { + bool isOperable() const { + if (!operable_under_only_bhp_limit) { + return false; + } else { + return ( (isOperableUnderBHPLimit() || isOperableUnderTHPLimit()) && + !(has_thp_constaint && !can_produce_inject_with_current_bhp) ); + } + } + + bool isOperableUnderBHPLimit() const { + return operable_under_only_bhp_limit && obey_thp_limit_under_bhp_limit; + } + + bool isOperableUnderTHPLimit() const { + return can_obtain_bhp_with_thp_limit && obey_bhp_limit_with_thp_limit; + } + + void reset() { + operable_under_only_bhp_limit = true; + obey_thp_limit_under_bhp_limit = true; + can_obtain_bhp_with_thp_limit = true; + obey_bhp_limit_with_thp_limit = true; + can_produce_inject_with_current_bhp = true; + has_thp_constaint = false; + } + + // whether the well can be operated under bhp limit + // without considering other limits. + // if it is false, then the well is not operable for sure. + bool operable_under_only_bhp_limit = true; + // if the well can be operated under bhp limit, will it obey(not violate) + // the thp limit when operated under bhp limit + bool obey_thp_limit_under_bhp_limit = true; + // whether the well operate under the thp limit only + bool can_obtain_bhp_with_thp_limit = true; + // whether the well obey bhp limit when operated under thp limit + bool obey_bhp_limit_with_thp_limit = true; + + // TODO: the following criterion is based on the current state of + // the well, we consider it is a numerical criterion. + // at the moment, we only apply it with well has THP constraint. + // whether the well can produce / inject with the current bhp of the well + // it might be updated with other criterion with investigation with more cases. + bool can_produce_inject_with_current_bhp = true; + // whether the well has a THP constraint + bool has_thp_constaint = false; + }; + } #include "WellInterface_impl.hpp" diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index cdb564ef3..80e81cb1d 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -457,8 +457,25 @@ namespace Opm if (wellhelpers::constraintBroken( well_state.bhp(), well_state.thp(), well_state.wellRates(), w, np, well_type_, wc, ctrl_index)) { - // ctrl_index will be the index of the broken constraint after the loop. - break; + + // if the well can not work under THP / BHP control, we should not switch to THP / BHP control + const bool cannot_switch_to_bhp = well_controls_iget_type(wc, ctrl_index) == BHP && !operability_status_.isOperableUnderBHPLimit(); + const bool cannot_switch_to_thp = well_controls_iget_type(wc, ctrl_index) == THP && !operability_status_.isOperableUnderTHPLimit(); + const bool cannot_switch = cannot_switch_to_bhp || cannot_switch_to_thp; + if ( !cannot_switch ) { + + // ctrl_index will be the index of the broken constraint after the loop. + break; + } else { + // before we figure out to handle it, we give some debug information here + if ( well_controls_iget_type(wc, ctrl_index) == BHP && !operability_status_.isOperableUnderBHPLimit() ) { + OpmLog::debug("well " + name() + " breaks the BHP limit, while it is not operable under BHP limit"); + } + + if ( well_controls_iget_type(wc, ctrl_index) == THP && !operability_status_.isOperableUnderTHPLimit() ) { + OpmLog::debug("well " + name() + " breaks the THP limit, while it is not operable under THP limit"); + } + } } } @@ -719,6 +736,15 @@ namespace Opm return; } + // Based on current understanding, only under prediction mode, we need to shut well due to various + // reasons or limits. With more knowlage or testing cases later, this might need to be corrected. + if (!underPredictionMode() ) { + return; + } + + // updating well test state based on physical (THP/BHP) limits. + updateWellTestStatePhysical(well_state, simulationTime, writeMessageToOPMLog, wellTestState); + // updating well test state based on Economic limits. updateWellTestStateEconomic(well_state, simulationTime, writeMessageToOPMLog, wellTestState); @@ -729,6 +755,30 @@ namespace Opm + template + void + WellInterface:: + updateWellTestStatePhysical(const WellState& well_state, + const double simulation_time, + const bool write_message_to_opmlog, + WellTestState& well_test_state) const + { + if (!isOperable()) { + well_test_state.addClosedWell(name(), WellTestConfig::Reason::PHYSICAL, simulation_time); + if (write_message_to_opmlog) { + // TODO: considering auto shut in? + const std::string msg = "well " + name() + + std::string(" will be shut as it can not operate under current reservoir condition"); + OpmLog::info(msg); + } + } + + } + + + + + template void WellInterface:: @@ -1207,4 +1257,15 @@ namespace Opm } + + + + template + bool + WellInterface:: + isOperable() const { + return operability_status_.isOperable(); + } + + }