From 94c0b49cf7b569f78f24e8bb257d552857248646 Mon Sep 17 00:00:00 2001 From: Stein Krogstad Date: Thu, 26 Oct 2023 17:28:05 +0200 Subject: [PATCH] revert changes --- opm/simulators/wells/MultisegmentWell.hpp | 9 +- .../wells/MultisegmentWell_impl.hpp | 101 +++++++- opm/simulators/wells/SingleWellState.cpp | 6 + opm/simulators/wells/SingleWellState.hpp | 4 + opm/simulators/wells/StandardWell.hpp | 9 +- opm/simulators/wells/StandardWell_impl.hpp | 227 ++++++++++++++--- opm/simulators/wells/VFPHelpers.cpp | 84 +++++++ opm/simulators/wells/VFPHelpers.hpp | 23 ++ opm/simulators/wells/VFPProdProperties.cpp | 19 +- opm/simulators/wells/VFPProdProperties.hpp | 6 +- opm/simulators/wells/VFPProperties.hpp | 17 ++ opm/simulators/wells/WellBhpThpCalculator.cpp | 116 +++++++++ opm/simulators/wells/WellBhpThpCalculator.hpp | 19 ++ opm/simulators/wells/WellInterface.hpp | 31 ++- opm/simulators/wells/WellInterfaceGeneric.cpp | 43 ++++ opm/simulators/wells/WellInterfaceGeneric.hpp | 5 + opm/simulators/wells/WellInterface_impl.hpp | 229 +++++++++++++++--- 17 files changed, 859 insertions(+), 89 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 09a019e3d..76c63303e 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -132,6 +132,8 @@ namespace Opm const WellState& well_state, DeferredLogger& deferred_logger) override; // should be const? + void updateIPRImplicit(const Simulator& ebos_simulator, WellState& well_state, DeferredLogger& deferred_logger); + virtual void updateProductivityIndex(const Simulator& ebosSimulator, const WellProdIndexCalculator& wellPICalc, WellState& well_state, @@ -246,6 +248,10 @@ namespace Opm const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const; + bool computeWellPotentialsImplicit(const Simulator& ebos_simulator, + std::vector& well_potentials, + DeferredLogger& deferred_logger) const; + virtual double getRefDensity() const override; virtual bool iterateWellEqWithControl(const Simulator& ebosSimulator, @@ -262,7 +268,8 @@ namespace Opm const Well::ProductionControls& prod_controls, WellState& well_state, const GroupState& group_state, - DeferredLogger& deferred_logger) override; + DeferredLogger& deferred_logger, + const bool allow_switch = true) override; virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator, const double dt, diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 7bdb507d5..325cd4639 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -286,13 +286,23 @@ namespace Opm } debug_cost_counter_ = 0; - // does the well have a THP related constraint? - const auto& summaryState = ebosSimulator.vanguard().summaryState(); - if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) { - computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger); - } else { - well_potentials = computeWellPotentialWithTHP( - well_state, ebosSimulator, deferred_logger); + bool converged_implicit = false; + if (this->param_.local_well_solver_control_switching_) { + converged_implicit = computeWellPotentialsImplicit(ebosSimulator, well_potentials, deferred_logger); + if (!converged_implicit) { + deferred_logger.debug("Implicit potential calculations failed for well " + + this->name() + ", reverting to original aproach."); + } + } + if (!converged_implicit) { + // does the well have a THP related constraint? + const auto& summaryState = ebosSimulator.vanguard().summaryState(); + if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) { + computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger); + } else { + well_potentials = computeWellPotentialWithTHP( + well_state, ebosSimulator, deferred_logger); + } } deferred_logger.debug("Cost in iterations of finding well potential for well " + this->name() + ": " + std::to_string(debug_cost_counter_)); @@ -488,7 +498,70 @@ namespace Opm return potentials; } + template + bool + MultisegmentWell:: + computeWellPotentialsImplicit(const Simulator& ebos_simulator, + std::vector& well_potentials, + DeferredLogger& deferred_logger) const + { + // Create a copy of the well. + // TODO: check if we can avoid taking multiple copies. Call from updateWellPotentials + // is allready a copy, but not from other calls. + MultisegmentWell well_copy(*this); + well_copy.debug_cost_counter_ = 0; + // store a copy of the well state, we don't want to update the real well state + WellState well_state_copy = ebos_simulator.problem().wellModel().wellState(); + const auto& group_state = ebos_simulator.problem().wellModel().groupState(); + auto& ws = well_state_copy.well(this->index_of_well_); + + // get current controls + const auto& summary_state = ebos_simulator.vanguard().summaryState(); + auto inj_controls = well_copy.well_ecl_.isInjector() + ? well_copy.well_ecl_.injectionControls(summary_state) + : Well::InjectionControls(0); + auto prod_controls = well_copy.well_ecl_.isProducer() + ? well_copy.well_ecl_.productionControls(summary_state) + : Well::ProductionControls(0); + + // prepare/modify well state and control + well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls); + + well_copy.scaleSegmentPressuresWithBhp(well_state_copy); + + // initialize rates from previous potentials + const int np = this->number_of_phases_; + bool trivial = true; + for (int phase = 0; phase < np; ++phase){ + trivial = trivial && (ws.well_potentials[phase] == 0.0) ; + } + if (!trivial) { + const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0; + for (int phase = 0; phase < np; ++phase) { + ws.surface_rates[phase] = sign * ws.well_potentials[phase]; + } + } + well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(), + this->segments_.perforations(), + well_state_copy); + + well_copy.calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger); + const double dt = ebos_simulator.timeStepSize(); + // solve equations + const bool converged = well_copy.iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, + deferred_logger); + + // fetch potentials (sign is updated on the outside). + well_potentials.clear(); + well_potentials.resize(np, 0.0); + for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) { + const EvalWell rate = well_copy.primary_variables_.getQs(compIdx); + well_potentials[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value(); + } + debug_cost_counter_ += well_copy.debug_cost_counter_; + return converged; + } template void @@ -1227,6 +1300,13 @@ namespace Opm } } + template + void + MultisegmentWell:: + updateIPRImplicit(const Simulator& ebos_simulator, WellState& well_state, DeferredLogger& deferred_logger) + { + } + template void MultisegmentWell:: @@ -1441,7 +1521,8 @@ namespace Opm const Well::ProductionControls& prod_controls, WellState& well_state, const GroupState& group_state, - DeferredLogger& deferred_logger) + DeferredLogger& deferred_logger, + const bool allow_switch /*true*/) { //if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true; @@ -1586,9 +1667,9 @@ namespace Opm } else { this->operability_status_.operable_under_only_bhp_limit = !is_stopped; } - // We reset the well status to it's original state. Status is updated + // We reset the well status to it's original state. Status is updated // on the outside based on operability status - this->wellStatus_ = well_status; + // this->wellStatus_ = well_status; } std::string message = fmt::format(" Well {} converged in {} inner iterations (" "{} control/status switches).", this->name(), it, switch_count); diff --git a/opm/simulators/wells/SingleWellState.cpp b/opm/simulators/wells/SingleWellState.cpp index 774fcdd39..975b7ee8d 100644 --- a/opm/simulators/wells/SingleWellState.cpp +++ b/opm/simulators/wells/SingleWellState.cpp @@ -42,6 +42,8 @@ SingleWellState::SingleWellState(const std::string& name_, , temperature(temp) , well_potentials(pu_.num_phases) , productivity_index(pu_.num_phases) + , ipr_a(pu_.num_phases) + , ipr_b(pu_.num_phases) , surface_rates(pu_.num_phases) , reservoir_rates(pu_.num_phases) , prev_surface_rates(pu_.num_phases) @@ -89,6 +91,8 @@ void SingleWellState::shut() { std::fill(this->prev_surface_rates.begin(), this->prev_surface_rates.end(), 0); std::fill(this->reservoir_rates.begin(), this->reservoir_rates.end(), 0); std::fill(this->productivity_index.begin(), this->productivity_index.end(), 0); + std::fill(this->ipr_a.begin(), this->ipr_a.end(), 0); + std::fill(this->ipr_b.begin(), this->ipr_b.end(), 0); auto& connpi = this->perf_data.prod_index; connpi.assign(connpi.size(), 0); @@ -305,6 +309,8 @@ bool SingleWellState::operator==(const SingleWellState& rhs) const this->phase_mixing_rates == rhs.phase_mixing_rates && this->well_potentials == rhs.well_potentials && this->productivity_index == rhs.productivity_index && + this->ipr_a == rhs.ipr_a && + this->ipr_b == rhs.ipr_b && this->surface_rates == rhs.surface_rates && this->reservoir_rates == rhs.reservoir_rates && this->prev_surface_rates == rhs.prev_surface_rates && diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp index b2e89585c..807df77c0 100644 --- a/opm/simulators/wells/SingleWellState.hpp +++ b/opm/simulators/wells/SingleWellState.hpp @@ -61,6 +61,8 @@ public: serializer(phase_mixing_rates); serializer(well_potentials); serializer(productivity_index); + serializer(ipr_a); + serializer(ipr_b); serializer(surface_rates); serializer(reservoir_rates); serializer(prev_surface_rates); @@ -98,6 +100,8 @@ public: std::vector well_potentials; std::vector productivity_index; + std::vector ipr_a; + std::vector ipr_b; std::vector surface_rates; std::vector reservoir_rates; std::vector prev_surface_rates; diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index a19379b7b..d5dc1a1f4 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -212,7 +212,8 @@ namespace Opm const Well::ProductionControls& prod_controls, WellState& well_state, const GroupState& group_state, - DeferredLogger& deferred_logger) override; + DeferredLogger& deferred_logger, + const bool allow_switch = true) override; /// \brief Wether the Jacobian will also have well contributions in it. virtual bool jacobianContainsWellContributions() const override @@ -240,6 +241,8 @@ namespace Opm const double alq_value, DeferredLogger& deferred_logger) const override; + void updateIPRImplicit(const Simulator& ebosSimulator, WellState& well_state, DeferredLogger& deferred_logger); + virtual void computeWellRatesWithBhp( const Simulator& ebosSimulator, const double& bhp, @@ -321,6 +324,10 @@ namespace Opm DeferredLogger& deferred_logger, const WellState &well_state) const; + bool computeWellPotentialsImplicit(const Simulator& ebos_simulator, + std::vector& well_potentials, + DeferredLogger& deferred_logger) const; + virtual double getRefDensity() const override; // get the mobility for specific perforation diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 03c5e980e..03f905bb5 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -832,6 +832,89 @@ namespace Opm this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size()); } + template + void + StandardWell:: + updateIPRImplicit(const Simulator& ebosSimulator, + WellState& well_state, + DeferredLogger& deferred_logger) + { + // Compute IPR based on *converged* well-equation: + // For a component rate r the derivative dr/dbhp is obtained by + // dr/dbhp = - (partial r/partial x) * inv(partial Eq/partial x) * (partial Eq/partial control_value) + // where Eq(x)=0 is the well equation setup with bhp control and primary varables x + //StandardWell well_copy(*this); + + // We shouldn't have zero rates at this stage, but check + bool zero_rates; + auto rates = well_state.well(this->index_of_well_).surface_rates; + zero_rates = true; + for (std::size_t p = 0; p < rates.size(); ++p) { + zero_rates &= rates[p] == 0.0; + } + auto& ws = well_state.well(this->index_of_well_); + if (zero_rates) { + const auto msg = fmt::format("updateIPRImplicit: Well {} has zero rate, reverting to explicit IPR-calulations", this->name()); + deferred_logger.debug(msg); + updateIPR(ebosSimulator, deferred_logger); + for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){ + const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx); + ws.ipr_a[idx] = this->ipr_a_[comp_idx]; + ws.ipr_b[idx] = this->ipr_b_[comp_idx]; + } + } else { + const auto& group_state = ebosSimulator.problem().wellModel().groupState(); + + // XXX maybe don't update this + std::fill(ws.ipr_a.begin(), ws.ipr_a.end(), 0.); + std::fill(ws.ipr_b.begin(), ws.ipr_b.end(), 0.); + //WellState well_state_copy = well_state; + auto inj_controls = Well::InjectionControls(0); + auto prod_controls = Well::ProductionControls(0); + prod_controls.addControl(Well::ProducerCMode::BHP); + prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp; + + // Set current control to bhp, and bhp value in state, modify bhp limit in control object. + const auto cmode = ws.production_cmode; + ws.production_cmode = Well::ProducerCMode::BHP; + const double dt = ebosSimulator.timeStepSize(); + assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + + const double nEq = this->primary_variables_.numWellEq(); + BVectorWell rhs(1); + rhs[0].resize(nEq); + // rhs = 0 except -1 for control eq + for (size_t i=0; i < nEq; ++i){ + rhs[0][i] = 0.0; + } + rhs[0][Bhp] = -1.0; + + BVectorWell x_well(1); + x_well[0].resize(nEq); + this->linSys_.solve(rhs, x_well); + + for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){ + EvalWell comp_rate = this->primary_variables_.getQs(comp_idx); + const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx); + for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) { + ws.ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq); + } + ws.ipr_a[idx] = ws.ipr_b[idx]*ws.bhp - comp_rate.value(); + + //for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) { + // this->ipr_b_[comp_idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq); + //} + // XXX maybe don't update this + //this->ipr_a_[comp_idx] = this->ipr_b_[comp_idx]*ws.bhp - comp_rate.value(); + // For ipr in well_state use same ordering as potentials etc. + //const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx); + //ws.ipr_a[idx] = this->ipr_a_[comp_idx]; + //ws.ipr_b[idx] = this->ipr_b_[comp_idx]; + } + // reset cmode + ws.production_cmode = cmode; + } + } template void @@ -1492,6 +1575,63 @@ namespace Opm return potentials; } + template + bool + StandardWell:: + computeWellPotentialsImplicit(const Simulator& ebos_simulator, + std::vector& well_potentials, + DeferredLogger& deferred_logger) const + { + // Create a copy of the well. + // TODO: check if we can avoid taking multiple copies. Call from updateWellPotentials + // is allready a copy, but not from other calls. + StandardWell well_copy(*this); + + // store a copy of the well state, we don't want to update the real well state + WellState well_state_copy = ebos_simulator.problem().wellModel().wellState(); + const auto& group_state = ebos_simulator.problem().wellModel().groupState(); + auto& ws = well_state_copy.well(this->index_of_well_); + + // get current controls + const auto& summary_state = ebos_simulator.vanguard().summaryState(); + auto inj_controls = well_copy.well_ecl_.isInjector() + ? well_copy.well_ecl_.injectionControls(summary_state) + : Well::InjectionControls(0); + auto prod_controls = well_copy.well_ecl_.isProducer() + ? well_copy.well_ecl_.productionControls(summary_state) : + Well::ProductionControls(0); + + // prepare/modify well state and control + well_copy.prepareForPotentialCalculations(summary_state, well_state_copy, inj_controls, prod_controls); + + // initialize rates from previous potentials + const int np = this->number_of_phases_; + bool trivial = true; + for (int phase = 0; phase < np; ++phase){ + trivial = trivial && (ws.well_potentials[phase] == 0.0) ; + } + if (!trivial) { + const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0; + for (int phase = 0; phase < np; ++phase) { + ws.surface_rates[phase] = sign * ws.well_potentials[phase]; + } + } + + well_copy.calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger); + const double dt = ebos_simulator.timeStepSize(); + // iterate to get a solution at the given bhp. + const bool converged = well_copy.iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, + deferred_logger); + + // fetch potentials (sign is updated on the outside). + well_potentials.clear(); + well_potentials.resize(np, 0.0); + for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) { + const EvalWell rate = well_copy.primary_variables_.getQs(compIdx); + well_potentials[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value(); + } + return converged; + } template @@ -1554,29 +1694,35 @@ namespace Opm return; } - // does the well have a THP related constraint? - const auto& summaryState = ebosSimulator.vanguard().summaryState(); - if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) { - // get the bhp value based on the bhp constraints - double bhp = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState); + bool converged_implicit = false; + if (this->param_.local_well_solver_control_switching_) { + converged_implicit = computeWellPotentialsImplicit(ebosSimulator, well_potentials, deferred_logger); + } + if (!converged_implicit) { + // does the well have a THP related constraint? + const auto& summaryState = ebosSimulator.vanguard().summaryState(); + if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) { + // get the bhp value based on the bhp constraints + double bhp = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState); - // In some very special cases the bhp pressure target are - // temporary violated. This may lead to too small or negative potentials - // that could lead to premature shutting of wells. - // As a remedy the bhp that gives the largest potential is used. - // For converged cases, ws.bhp <=bhp for injectors and ws.bhp >= bhp, - // and the potentials will be computed using the limit as expected. - const auto& ws = well_state.well(this->index_of_well_); - if (this->isInjector()) - bhp = std::max(ws.bhp, bhp); - else - bhp = std::min(ws.bhp, bhp); + // In some very special cases the bhp pressure target are + // temporary violated. This may lead to too small or negative potentials + // that could lead to premature shutting of wells. + // As a remedy the bhp that gives the largest potential is used. + // For converged cases, ws.bhp <=bhp for injectors and ws.bhp >= bhp, + // and the potentials will be computed using the limit as expected. + const auto& ws = well_state.well(this->index_of_well_); + if (this->isInjector()) + bhp = std::max(ws.bhp, bhp); + else + bhp = std::min(ws.bhp, bhp); - assert(std::abs(bhp) != std::numeric_limits::max()); - computeWellRatesWithBhpIterations(ebosSimulator, bhp, well_potentials, deferred_logger); - } else { - // the well has a THP related constraint - well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger, well_state); + assert(std::abs(bhp) != std::numeric_limits::max()); + computeWellRatesWithBhpIterations(ebosSimulator, bhp, well_potentials, deferred_logger); + } else { + // the well has a THP related constraint + well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger, well_state); + } } this->checkNegativeWellPotentials(well_potentials, @@ -2174,32 +2320,40 @@ namespace Opm const Well::ProductionControls& prod_controls, WellState& well_state, const GroupState& group_state, - DeferredLogger& deferred_logger) + DeferredLogger& deferred_logger, + const bool allow_switch /*true*/) { const int max_iter = this->param_.max_inner_iter_wells_; - + int it = 0; bool converged; bool relax_convergence = false; this->regularize_ = false; const auto& summary_state = ebosSimulator.vanguard().summaryState(); - // Max status switch frequency should be 2 to avoid getting stuck in cycle - constexpr int min_its_after_switch = 2; + // Max status switch frequency should be 2 to avoid getting stuck in cycle + constexpr int min_its_after_switch = 4; int its_since_last_switch = min_its_after_switch; int switch_count= 0; - const auto well_status = this->wellStatus_; + // if we fail to solve eqs, we reset status/control before leaving + const auto well_status_orig = this->wellStatus_; + auto well_status_cur = well_status_orig; + int status_switch_count = 0; const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN); bool changed = false; - bool final_check = false; + bool final_check = false; do { its_since_last_switch++; if (allow_switching && its_since_last_switch >= min_its_after_switch){ const double wqTotal = this->primary_variables_.eval(WQTotal).value(); - changed = this->updateWellControlAndStatusLocalIteration(ebosSimulator, well_state, group_state, inj_controls, prod_controls, wqTotal, deferred_logger); + changed = this->updateWellControlAndStatusLocalIteration(ebosSimulator, well_state, group_state, inj_controls, prod_controls, wqTotal, deferred_logger, allow_switch); if (changed){ its_since_last_switch = 0; switch_count++; + if (well_status_cur != this->wellStatus_) { + well_status_cur = this->wellStatus_; + status_switch_count++; + } } if (!changed && final_check) { break; @@ -2207,7 +2361,7 @@ namespace Opm final_check = false; } } - + assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); if (it > this->param_.strict_inner_iter_wells_) { @@ -2226,14 +2380,20 @@ namespace Opm its_since_last_switch = min_its_after_switch; } else { break; - } + } } ++it; solveEqAndUpdateWellState(summary_state, well_state, deferred_logger); initPrimaryVariablesEvaluation(); - } while (it < max_iter); + if (status_switch_count > 10) { + // constantly changing status, give up (do operability check on the outside) + converged = false; + break; + } + } while (it < max_iter); + if (converged) { if (allow_switching){ // update operability if status change @@ -2251,11 +2411,12 @@ namespace Opm // For this to happen, isOperableAndSolvable() must change from true to false, // and (until the most recent commit) the well needs to be open for this to trigger. // Hence, the resetting of status. - this->wellStatus_ = well_status; + //this->wellStatus_ = well_status; } } else { + this->wellStatus_ = well_status_orig; const std::string message = fmt::format(" Well {} did not converged in {} inner iterations (" - "{} control/status switches).", this->name(), it, switch_count); + "{} switches, {} status changes).", this->name(), it, switch_count, status_switch_count); deferred_logger.debug(message); // add operability here as well ? } diff --git a/opm/simulators/wells/VFPHelpers.cpp b/opm/simulators/wells/VFPHelpers.cpp index 5ce6c42f4..1c6074543 100644 --- a/opm/simulators/wells/VFPHelpers.cpp +++ b/opm/simulators/wells/VFPHelpers.cpp @@ -499,6 +499,90 @@ double findTHP(const std::vector& bhp_array, return thp; } +std::pair +getMinimumBHPCoordinate(const VFPProdTable& table, + const double thp, + const double wfr, + const double gfr, + const double alq) +{ + double flo_at_bhp_min = 0.0; // start by checking flo=0 + auto flo_i = detail::findInterpData(flo_at_bhp_min, table.getFloAxis()); + auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); + auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis()); + auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis()); + auto alq_i = detail::findInterpData( alq, table.getALQAxis()); + + detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + double bhp_min = bhp_i.value; + std::vector flos = table.getFloAxis(); + for (size_t i = 0; i < flos.size(); ++i) { + flo_i = detail::findInterpData(flos[i], table.getFloAxis()); + bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + if (bhp_i.value < bhp_min){ + bhp_min = bhp_i.value; + flo_at_bhp_min = flos[i]; + } + } + // return negative flo + return std::make_pair(-flo_at_bhp_min, bhp_min); +} + +std::optional> +intersectWithIPR(const VFPProdTable& table, + const double thp, + const double wfr, + const double gfr, + const double alq, + const double ipr_a, + const double ipr_b) +{ + // NOTE: ipr-line is q=b*bhp - a! + // ipr is given for negative flo, so + // flo = -b*bhp + a, i.e., bhp = -(flo-a)/b + auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); + auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis()); + auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis()); + auto alq_i = detail::findInterpData( alq, table.getALQAxis()); + + if (ipr_b == 0.0) { + // this shouldn't happen, but deal with it to be safe + auto flo_i = detail::findInterpData(-ipr_a, table.getFloAxis()); + detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + return std::make_pair(ipr_b, bhp_i.value); + } + // find largest flo (flo_x) for which y = bhp(flo) + (flo-a)/b = 0 and dy/dflo > 0 + double flo_x = -1.0; + double flo0, flo1; + double y0, y1; + flo0 = 0.0; // start by checking flo=0 + auto flo_i = detail::findInterpData(flo0, table.getFloAxis()); + detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + y0 = bhp_i.value - ipr_a/ipr_b; // +0.0/ipr_b + + std::vector flos = table.getFloAxis(); + for (size_t i = 0; i < flos.size(); ++i) { + flo1 = flos[i]; + flo_i = detail::findInterpData(flo1, table.getFloAxis()); + bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + y1 = bhp_i.value + (flo1 - ipr_a)/ipr_b; + if (y0 < 0 && y1 >= 0){ + // crossing with positive slope + double w = -y0/(y1-y0); + w = std::clamp(w, 0.0, 1.0); // just to be safe (if y0~y1) + flo_x = flo0 + w*(flo1 - flo0); + } + flo0 = flo1; + y0 = y1; + } + // return (last) intersection if found (negative flo) + if (flo_x >= 0) { + return std::make_pair(-flo_x, -(flo_x - ipr_a)/ipr_b); + } else { + return std::nullopt; + } +} + template T getFlo(const VFPProdTable& table, const T& aqua, diff --git a/opm/simulators/wells/VFPHelpers.hpp b/opm/simulators/wells/VFPHelpers.hpp index 31b4702a3..b69bc1c22 100644 --- a/opm/simulators/wells/VFPHelpers.hpp +++ b/opm/simulators/wells/VFPHelpers.hpp @@ -25,6 +25,7 @@ #include #include #include +#include /** * This file contains a set of helper functions used by VFPProd / VFPInj. @@ -181,6 +182,28 @@ double findTHP(const std::vector& bhp_array, const std::vector& thp_array, double bhp); +/** +* Get (flo, bhp) at minimum bhp for given thp,wfr,gfr,alq +*/ +std::pair +getMinimumBHPCoordinate(const VFPProdTable& table, + const double thp, + const double wfr, + const double gfr, + const double alq); + +/** +* Get (flo, bhp) at largest occuring stable vfp/ipr-intersection +* if it exists +*/ +std::optional> +intersectWithIPR(const VFPProdTable& table, + const double thp, + const double wfr, + const double gfr, + const double alq, + const double ipr_a, + const double ipr_b); } // namespace detail diff --git a/opm/simulators/wells/VFPProdProperties.cpp b/opm/simulators/wells/VFPProdProperties.cpp index 5c9f9f420..f69df3572 100644 --- a/opm/simulators/wells/VFPProdProperties.cpp +++ b/opm/simulators/wells/VFPProdProperties.cpp @@ -137,6 +137,22 @@ bhpwithflo(const std::vector& flos, return bhps; } +double +VFPProdProperties:: +minimumBHP(const int table_id, + const double thp, + const double wfr, + const double gfr, + const double alq) const +{ + // Get the table + const VFPProdTable& table = detail::getTable(m_tables, table_id); + const auto retval = detail::getMinimumBHPCoordinate(table, thp, wfr, gfr, alq); + // returned pair is (flo, bhp) + return retval.second; +} + + void VFPProdProperties::addTable(const VFPProdTable& new_table) { this->m_tables.emplace( new_table.getTableNum(), new_table ); @@ -176,7 +192,8 @@ EvalWell VFPProdProperties::bhp(const int table_id, detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); - bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (std::max(0.0, bhp_val.dflo) * flo); + //bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (std::max(0.0, bhp_val.dflo) * flo); + bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (bhp_val.dflo* flo); bhp.setValue(bhp_val.value); return bhp; diff --git a/opm/simulators/wells/VFPProdProperties.hpp b/opm/simulators/wells/VFPProdProperties.hpp index 4dbb8ca7d..d5a72975a 100644 --- a/opm/simulators/wells/VFPProdProperties.hpp +++ b/opm/simulators/wells/VFPProdProperties.hpp @@ -129,7 +129,11 @@ public: return m_tables.empty(); } - + /** + * Returns minimum bhp for given thp, wfr, gfr and alq + */ + double minimumBHP(const int table_id, const double thp, + const double wfr, const double gfr, const double alq) const; protected: // calculate a group bhp values with a group of flo rate values std::vector bhpwithflo(const std::vector& flos, diff --git a/opm/simulators/wells/VFPProperties.hpp b/opm/simulators/wells/VFPProperties.hpp index 2f41d4881..8d86a4897 100644 --- a/opm/simulators/wells/VFPProperties.hpp +++ b/opm/simulators/wells/VFPProperties.hpp @@ -92,6 +92,23 @@ public: return detail::getGFR(table, aqua, liquid, vapour); } + std::pair + getFloIPR(const int table_id, const std::size_t well_index) const { + std::pairretval(0.0, 0.0); + const VFPProdTable& table = this->m_prod.getTable(table_id); + const auto& pu = well_state_.phaseUsage(); + const auto& ipr_a= well_state_.well(well_index).ipr_a; + const auto& aqua_a = pu.phase_used[BlackoilPhases::Aqua]? ipr_a[pu.phase_pos[BlackoilPhases::Aqua]]:0.0; + const auto& liquid_a = pu.phase_used[BlackoilPhases::Liquid]? ipr_a[pu.phase_pos[BlackoilPhases::Liquid]]:0.0; + const auto& vapour_a = pu.phase_used[BlackoilPhases::Vapour]? ipr_a[pu.phase_pos[BlackoilPhases::Vapour]]:0.0; + const auto& ipr_b= well_state_.well(well_index).ipr_b; + const auto& aqua_b = pu.phase_used[BlackoilPhases::Aqua]? ipr_b[pu.phase_pos[BlackoilPhases::Aqua]]:0.0; + const auto& liquid_b = pu.phase_used[BlackoilPhases::Liquid]? ipr_b[pu.phase_pos[BlackoilPhases::Liquid]]:0.0; + const auto& vapour_b = pu.phase_used[BlackoilPhases::Vapour]? ipr_b[pu.phase_pos[BlackoilPhases::Vapour]]:0.0; + return std::make_pair(detail::getFlo(table, aqua_a, liquid_a, vapour_a), + detail::getFlo(table, aqua_b, liquid_b, vapour_b)); + } + private: VFPInjProperties m_inj; VFPProdProperties m_prod; diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index 52ad20be2..b54d679c5 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -387,6 +387,32 @@ calculateBhpFromThp(const WellState& well_state, return bhp_tab - dp_hydro + bhp_adjustment; } +double +WellBhpThpCalculator:: +calculateMinimumBhpFromThp(const WellState& well_state, + const Well& well, + const SummaryState& summaryState, + const double rho) const +{ + assert(well_.isProducer()); // only producers can go here for now + + const double thp_limit = well_.getTHPConstraint(summaryState); + + const auto& controls = well.productionControls(summaryState); + const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell()); + const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell()); + + const double bhp_min = well_.vfpProperties()->getProd()->minimumBHP(controls.vfp_table_number, + thp_limit, wfr, gfr, + well_.getALQ(well_state)); + + const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth(); + const auto bhp_adjustment = getVfpBhpAdjustment(bhp_min, thp_limit); + const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, + rho, well_.gravity()); + return bhp_min - dp_hydro + bhp_adjustment; +} + double WellBhpThpCalculator::getVfpBhpAdjustment(const double bhp_tab, const double thp_limit) const { return well_.wellEcl().getWVFPDP().getPressureLoss(bhp_tab, thp_limit); @@ -839,6 +865,96 @@ bruteForceBracket(const std::function& eq, return bracket_found; } +bool WellBhpThpCalculator:: +isStableSolution(const WellState& well_state, + const Well& well, + const std::vector& rates, + const SummaryState& summaryState, + DeferredLogger& deferred_logger) const +{ + assert(int(rates.size()) == 3); // the vfp related only supports three phases now. + assert(well_.isProducer()); // only valid for producers + + static constexpr int Gas = BlackoilPhases::Vapour; + static constexpr int Oil = BlackoilPhases::Liquid; + static constexpr int Water = BlackoilPhases::Aqua; + + const double aqua = rates[Water]; + const double liquid = rates[Oil]; + const double vapour = rates[Gas]; + const double thp = well_.getTHPConstraint(summaryState); + + const auto& controls = well.productionControls(summaryState); + const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell()); + const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell()); + + const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number); + const bool use_vfpexplicit = well_.useVfpExplicit(); + + const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth(); + const auto bhp_adjustment = getVfpBhpAdjustment(well_state.well(well_.indexOfWell()).bhp, thp); + // XXX this needs to be fixed + //assert(bhp_adjustment == 0.0); + + const detail::VFPEvaluation bhp = detail::bhp(table, aqua, liquid, vapour, thp, well_.getALQ(well_state), wfr, gfr, use_vfpexplicit); + + if (bhp.dflo >= 0) { + return true; + } else { // maybe check if ipr is available + const auto ipr = well_.vfpProperties()->getFloIPR(controls.vfp_table_number, well_.indexOfWell()); + return bhp.dflo + 1/ipr.second >= 0; + } +} + +std::optional WellBhpThpCalculator:: +estimateStableBhp(const WellState& well_state, + const Well& well, + const std::vector& rates, + const double rho, + const SummaryState& summaryState, + DeferredLogger& deferred_logger) const +{ + // Given a *converged* well_state with ipr, estimate bhp of the stable solution + const auto& controls = well.productionControls(summaryState); + const double thp = well_.getTHPConstraint(summaryState); + const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number); + + const double aqua = rates[BlackoilPhases::Aqua]; + const double liquid = rates[BlackoilPhases::Liquid]; + const double vapour = rates[BlackoilPhases::Vapour]; + double flo = detail::getFlo(table, aqua, liquid, vapour); + double wfr, gfr; + if (well_.useVfpExplicit() || -flo < table.getFloAxis().front()) { + wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell()); + gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell()); + } else { + wfr = detail::getWFR(table, aqua, liquid, vapour); + gfr = detail::getGFR(table, aqua, liquid, vapour); + } + if (wfr <= 0.0 && gfr <= 0.0) { + // warning message + } + auto ipr = well_.vfpProperties()->getFloIPR(controls.vfp_table_number, well_.indexOfWell()); + + const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth(); + const auto bhp_adjustment = getVfpBhpAdjustment(well_state.well(well_.indexOfWell()).bhp, thp); + // XXX this needs to be fixed + //assert(bhp_adjustment == 0.0); + + const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, + rho, well_.gravity()); + if (ipr.first <= 0.0) { + // error message + } + const auto retval = detail::intersectWithIPR(table, thp, wfr, gfr, well_.getALQ(well_state), ipr.first+ipr.second*dp_hydro, ipr.second); + if (retval.has_value()) { + // returned pair is (flo, bhp) + return retval.value().second - dp_hydro; + } else { + return std::nullopt; + } +} + #define INSTANCE(...) \ template __VA_ARGS__ WellBhpThpCalculator:: \ calculateBhpFromThp<__VA_ARGS__>(const WellState&, \ diff --git a/opm/simulators/wells/WellBhpThpCalculator.hpp b/opm/simulators/wells/WellBhpThpCalculator.hpp index d762dd2c7..628e73a7b 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.hpp +++ b/opm/simulators/wells/WellBhpThpCalculator.hpp @@ -98,6 +98,25 @@ public: const double rho, DeferredLogger& deferred_logger) const; + double calculateMinimumBhpFromThp(const WellState& well_state, + const Well& well, + const SummaryState& summaryState, + const double rho) const; + + bool isStableSolution(const WellState& well_state, + const Well& well, + const std::vector& rates, + const SummaryState& summaryState, + DeferredLogger& deferred_logger) const; + + std::optional + estimateStableBhp (const WellState& well_state, + const Well& well, + const std::vector& rates, + const double rho, + const SummaryState& summaryState, + DeferredLogger& deferred_logger) const; + private: //! \brief Compute BHP from THP limit for an injector - implementation. template diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 3f914ea3f..7a24811c0 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -247,7 +247,8 @@ public: const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, const double WQTotal, - DeferredLogger& deferred_logger); + DeferredLogger& deferred_logger, + const bool allow_switch=true); virtual void updatePrimaryVariables(const SummaryState& summary_state, const WellState& well_state, @@ -414,7 +415,12 @@ protected: const WellProductionControls& prod_controls, WellState& well_state, const GroupState& group_state, - DeferredLogger& deferred_logger) = 0; + DeferredLogger& deferred_logger, + const bool allow_switch = true) = 0; + + virtual void updateIPRImplicit(const Simulator& ebosSimulator, + WellState& well_state, + DeferredLogger& deferred_logger) = 0; bool iterateWellEquations(const Simulator& ebosSimulator, const double dt, @@ -422,6 +428,27 @@ protected: const GroupState& group_state, DeferredLogger& deferred_logger); + bool robustWellSolveWithTHP(const Simulator& ebos_simulator, + const double dt, + const Well::InjectionControls& inj_controls, + const Well::ProductionControls& prod_controls, + WellState& well_state, + const GroupState& group_state, + DeferredLogger& deferred_logger); + + std::optional estimateOperableBhp(const Simulator& ebos_simulator, + const double dt, + WellState& well_state, + const GroupState& group_state, + const SummaryState& summary_state, + DeferredLogger& deferred_logger); + + bool solveWellWithBhp(const Simulator& ebos_simulator, + const double dt, + const double bhp, + WellState& well_state, + DeferredLogger& deferred_logger); + bool solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const GroupState& group_state, DeferredLogger& deferred_logger); diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 861bbbaa9..1f4676c0d 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -744,4 +744,47 @@ checkNegativeWellPotentials(std::vector& well_potentials, } } +void WellInterfaceGeneric:: +prepareForPotentialCalculations(const SummaryState& summary_state, + WellState& well_state, + Well::InjectionControls& inj_controls, + Well::ProductionControls& prod_controls) +{ + const bool has_thp = this->wellHasTHPConstraints(summary_state); + auto& ws = well_state.well(this->index_of_well_); + // Modify control (only pressure constraints) and set new target if needed. + // Also set value for current target in state + if (this->isInjector()) { + inj_controls.clearControls(); + inj_controls.addControl(Well::InjectorCMode::BHP); + if (has_thp){ + inj_controls.addControl(Well::InjectorCMode::THP); + } + if (!(ws.injection_cmode == Well::InjectorCMode::BHP)){ + if (has_thp){ + ws.injection_cmode = Well::InjectorCMode::THP; + ws.thp = this->getTHPConstraint(summary_state); + } else { + ws.injection_cmode = Well::InjectorCMode::BHP; + ws.bhp = inj_controls.bhp_limit; + } + } + } else { // producer + prod_controls.clearControls(); + prod_controls.addControl(Well::ProducerCMode::BHP); + if (has_thp){ + prod_controls.addControl(Well::ProducerCMode::THP); + } + if (!(ws.production_cmode == Well::ProducerCMode::BHP)){ + if (has_thp){ + ws.production_cmode = Well::ProducerCMode::THP; + ws.thp = this->getTHPConstraint(summary_state); + } else { + ws.production_cmode = Well::ProducerCMode::BHP; + ws.bhp = prod_controls.bhp_limit; + } + } + } +} + } // namespace Opm diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 0ec167046..966336989 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -245,6 +245,11 @@ protected: const bool checkOperability, DeferredLogger& deferred_logger); + void prepareForPotentialCalculations(const SummaryState& summary_state, + WellState& well_state, + Well::InjectionControls& inj_controls, + Well::ProductionControls& prod_controls); + // definition of the struct OperabilityStatus struct OperabilityStatus { bool isOperableAndSolvable() const { diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 4f7f8cfdc..c5a97092f 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -264,11 +264,12 @@ namespace Opm const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, const double wqTotal, - DeferredLogger& deferred_logger) + DeferredLogger& deferred_logger, + const bool allow_switch /*true*/) { const auto& summary_state = ebos_simulator.vanguard().summaryState(); const auto& schedule = ebos_simulator.vanguard().schedule(); - + if (this->wellUnderZeroRateTarget(summary_state, well_state) || !(this->well_ecl_.getStatus() == WellStatus::OPEN)) { return false; } @@ -276,30 +277,33 @@ namespace Opm const double sgn = this->isInjector() ? 1.0 : -1.0; if (!this->wellIsStopped()){ if (wqTotal*sgn <= 0.0){ + //std::cout << "Stopping well:" << this->name() << std::endl; this->stopWell(); return true; } else { bool changed = false; - auto& ws = well_state.well(this->index_of_well_); - const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) : - prod_controls.hasControl(Well::ProducerCMode::GRUP); + if (allow_switch) { + auto& ws = well_state.well(this->index_of_well_); + const bool hasGroupControl = this->isInjector() ? inj_controls.hasControl(Well::InjectorCMode::GRUP) : + prod_controls.hasControl(Well::ProducerCMode::GRUP); - changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls); - // TODO: with current way, the checkGroupConstraints might overwrite the result from checkIndividualConstraints, which remains to be investigated - if (hasGroupControl) { - changed = this->checkGroupConstraints(well_state, group_state, schedule, summary_state,deferred_logger); - } - - if (changed) { - const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP : - ws.production_cmode == Well::ProducerCMode::THP; - if (!thp_controlled){ - // don't call for thp since this might trigger additional local solve - updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger); - } else { - ws.thp = this->getTHPConstraint(summary_state); + changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls); + // TODO: with current way, the checkGroupConstraints might overwrite the result from checkIndividualConstraints, which remains to be investigated + if (hasGroupControl) { + changed = changed || this->checkGroupConstraints(well_state, group_state, schedule, summary_state,deferred_logger); + } + + if (changed) { + const bool thp_controlled = this->isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP : + ws.production_cmode == Well::ProducerCMode::THP; + if (!thp_controlled){ + // don't call for thp since this might trigger additional local solve + updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger); + } else { + ws.thp = this->getTHPConstraint(summary_state); + } + updatePrimaryVariables(summary_state, well_state, deferred_logger); } - updatePrimaryVariables(summary_state, well_state, deferred_logger); } return changed; } @@ -311,20 +315,29 @@ namespace Opm const bool has_thp = this->wellHasTHPConstraints(summary_state); if (has_thp){ // calculate bhp from thp-limit (using explicit fractions zince zero rate) - // TODO: this will often be too strict condition for re-opening, a better - // option is probably minimum bhp on current vfp-curve, but some more functionality - // is needed for this option to be robustly implemented. + // TODO: this will often be too strict condition for re-opening, a better + // option is probably minimum bhp on current vfp-curve, but some more functionality + // is needed for this option to be robustly implemented. std::vector rates(this->num_components_); - const double bhp_thp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state, rates, this->well_ecl_, summary_state, this->getRefDensity(), deferred_logger); + //const double bhp_thp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state, rates, this->well_ecl_, summary_state, this->getRefDensity(), deferred_logger); if (this->isInjector()){ + const double bhp_thp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state, rates, this->well_ecl_, summary_state, this->getRefDensity(), deferred_logger); inj_limit = std::min(bhp_thp, inj_controls.bhp_limit); } else { - prod_limit = std::max(bhp_thp, prod_controls.bhp_limit); + const double bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity()); + prod_limit = std::max(bhp_min, prod_controls.bhp_limit); + //auto prates = well_state.well(this->index_of_well_).prev_surface_rates; + //std::cout << this->name() << ": Min bhp: " << bhp_min << " prod limit: " << prod_limit << " prev rates: " << prates[0] << " " << prates[1] << " " << prates[2] << std::endl; } } const double bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit; if (bhp_diff > 0){ + //std::cout << "Re-opening well:" << this->name() << std::endl; this->openWell(); + well_state.well(this->index_of_well_).bhp = prod_limit; + if (has_thp) { + well_state.well(this->index_of_well_).thp = this->getTHPConstraint(summary_state); + } return true; } else { return false; @@ -346,7 +359,9 @@ namespace Opm WellState well_state_copy = well_state; auto& ws = well_state_copy.well(this->indexOfWell()); - + if (ws.production_cmode == Well::ProducerCMode::GRUP) { + ws.production_cmode = Well::ProducerCMode::BHP; + } updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger); calculateExplicitQuantities(simulator, well_state_copy, deferred_logger); const auto& summary_state = simulator.vanguard().summaryState(); @@ -354,12 +369,7 @@ namespace Opm initPrimaryVariablesEvaluation(); if (this->isProducer()) { - const auto& schedule = simulator.vanguard().schedule(); - const auto report_step = simulator.episodeIndex(); - const auto& glo = schedule.glo(report_step); - if (glo.active()) { - gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger); - } + gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger); } WellTestState welltest_state_temp; @@ -448,9 +458,13 @@ namespace Opm if (!this->param_.local_well_solver_control_switching_){ converged = this->iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); } else { - converged = this->iterateWellEqWithSwitching(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) { + converged = robustWellSolveWithTHP(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + } else { + converged = this->iterateWellEqWithSwitching(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + } } - + } catch (NumericalProblem& e ) { const std::string msg = "Inner well iterations failed for well " + this->name() + " Treat the well as unconverged. "; deferred_logger.warning("INNER_ITERATION_FAILED", msg); @@ -459,6 +473,141 @@ namespace Opm return converged; } + template + bool + WellInterface:: + robustWellSolveWithTHP(const Simulator& ebos_simulator, + const double dt, + const Well::InjectionControls& inj_controls, + const Well::ProductionControls& prod_controls, + WellState& well_state, + const GroupState& group_state, + DeferredLogger& deferred_logger) + { + const auto& summary_state = ebos_simulator.vanguard().summaryState(); + bool is_operable = true; + bool converged; + converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + if (converged && !this->stopppedOrZeroRateTarget(summary_state, well_state) && well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP) { + auto rates = well_state.well(this->index_of_well_).surface_rates; + this->adaptRatesForVFP(rates); + // FIX !!!!!!!!!! + // first, check just dbhp/dflo + + // update ipr for stability-check + this->updateIPRImplicit(ebos_simulator, well_state, deferred_logger); + bool is_stable = WellBhpThpCalculator(*this).isStableSolution(well_state, this->well_ecl_, rates, summary_state, deferred_logger); + if (!is_stable) { + // solution converged to an unstable point! + //this->operability_status_.use_vfpexplicit = true; + auto bhp_stable = WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state, deferred_logger); + if (!bhp_stable.has_value()) { + // well is not operable (no intersection) + is_operable = false; + } else { + // solve equations for bhp_stable (in attempt to get closer to actual solution) + const bool converged_bhp = solveWellWithBhp(ebos_simulator, dt, bhp_stable.value(), well_state, deferred_logger); + if (!converged_bhp) { + // + } else { + // re-solve well equations + converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + } + } + } + } else if (!converged) { + // Well did not converge, switch to explicit fractions + this->operability_status_.use_vfpexplicit = true; + //auto well_state_copy = well_state; + // + auto bhp_target = estimateOperableBhp(ebos_simulator, dt, well_state, group_state, summary_state, deferred_logger); + if (!bhp_target.has_value()) { + // well can't operate using explicit fractions + is_operable = false; + } else { + // first solve well equations for bhp = bhp_target which should give solution close to actual + converged = solveWellWithBhp(ebos_simulator, dt, bhp_target.value(), well_state, deferred_logger); + if (!converged) { + // debug + } else { + // re-solve well equations + converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + } + } + } + if (!is_operable) { + this->stopWell(); + if (this->wellHasTHPConstraints(summary_state)){ + this->operability_status_.can_obtain_bhp_with_thp_limit = false; + this->operability_status_.obey_thp_limit_under_bhp_limit = false; + } else { + this->operability_status_.operable_under_only_bhp_limit = false; + } + } + return converged; + } + + template + std::optional + WellInterface:: + estimateOperableBhp(const Simulator& ebos_simulator, + const double dt, + WellState& well_state, + const GroupState& group_state, + const SummaryState& summary_state, + DeferredLogger& deferred_logger) + { + // Given an unconverged well or closed well, estimate an operable bhp (if any) + // Get minimal bhp from vfp-curve + double bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity()); + // Solve + const bool converged = solveWellWithBhp(ebos_simulator, dt, bhp_min, well_state, deferred_logger); + if (!converged) { + //report problem - return null + } else if (this->wellIsStopped()) { + return std::nullopt; + } + this->updateIPRImplicit(ebos_simulator, well_state, deferred_logger); + auto rates = well_state.well(this->index_of_well_).surface_rates; + this->adaptRatesForVFP(rates); + return WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state, deferred_logger); + } + + template + bool + WellInterface:: + solveWellWithBhp(const Simulator& ebos_simulator, + const double dt, + const double bhp, + WellState& well_state, + DeferredLogger& deferred_logger) + { + // Solve a well using single bhp-constraint (but close if not operable under this) + auto group_state = GroupState(); // empty group + auto inj_controls = Well::InjectionControls(0); + auto prod_controls = Well::ProductionControls(0); + auto& ws = well_state.well(this->index_of_well_); + auto cmode_inj = ws.injection_cmode; + auto cmode_prod = ws.production_cmode; + if (this->isInjector()) { + inj_controls.addControl(Well::InjectorCMode::BHP); + inj_controls.bhp_limit = bhp; + inj_controls.cmode = Well::InjectorCMode::BHP; + ws.injection_cmode = Well::InjectorCMode::BHP; + } else { + prod_controls.addControl(Well::ProducerCMode::BHP); + prod_controls.bhp_limit = bhp; + prod_controls.cmode = Well::ProducerCMode::BHP; + ws.production_cmode = Well::ProducerCMode::BHP; + } + // update well-state + ws.bhp = bhp; + // solve + const bool converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger, /*allow_switch*/false); + ws.injection_cmode = cmode_inj; + ws.production_cmode = cmode_prod; + return converged; + } template bool @@ -717,8 +866,8 @@ namespace Opm const auto& schedule = ebos_simulator.vanguard().schedule(); auto report_step_idx = ebos_simulator.episodeIndex(); const auto& glo = schedule.glo(report_step_idx); - if(glo.active() && glo.has_well(well_name)) { - const auto increment = glo.gaslift_increment(); + if(glo.has_well(well_name)) { + auto increment = glo.gaslift_increment(); auto alq = well_state.getALQ(well_name); bool converged; while (alq > 0) { @@ -751,8 +900,9 @@ namespace Opm deferred_logger.info(msg); return; } + const auto& well_ecl = this->wellEcl(); const auto& schedule = ebos_simulator.vanguard().schedule(); - const auto report_step_idx = ebos_simulator.episodeIndex(); + auto report_step_idx = ebos_simulator.episodeIndex(); const auto& glo = schedule.glo(report_step_idx); if (!glo.has_well(well_name)) { const std::string msg = fmt::format( @@ -768,7 +918,6 @@ namespace Opm max_alq = *max_alq_optional; } else { - const auto& well_ecl = this->wellEcl(); const auto& controls = well_ecl.productionControls(summary_state); const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number); const auto& alq_values = table.getALQAxis(); @@ -787,7 +936,7 @@ namespace Opm updateWellOperability(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger) - { + { if (this->param_.local_well_solver_control_switching_) { const bool success = updateWellOperabilityFromWellEq(ebos_simulator, well_state, deferred_logger); if (success) { @@ -833,7 +982,7 @@ namespace Opm // equations should be converged at this stage, so only one it is needed bool converged = iterateWellEquations(ebos_simulator, dt, well_state_copy, group_state, deferred_logger); return converged; - } + } template void