From 7c91c015cfc1f9ba6c9410539db879f035eb7e56 Mon Sep 17 00:00:00 2001 From: Stein Krogstad Date: Thu, 7 Dec 2023 12:15:40 +0100 Subject: [PATCH] updates according to Atgeirrs comments --- .../wells/MultisegmentWell_impl.hpp | 11 ++- opm/simulators/wells/StandardWell.hpp | 2 +- opm/simulators/wells/StandardWell_impl.hpp | 22 +++--- opm/simulators/wells/VFPHelpers.cpp | 2 +- opm/simulators/wells/WellBhpThpCalculator.cpp | 2 +- opm/simulators/wells/WellInterfaceGeneric.cpp | 2 +- opm/simulators/wells/WellInterfaceGeneric.hpp | 2 +- opm/simulators/wells/WellInterface_impl.hpp | 78 +++++++++---------- 8 files changed, 61 insertions(+), 60 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index d213d96d1..eda6ae209 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1312,8 +1312,8 @@ namespace Opm { // 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 + // dr/dbhp = - (partial r/partial x) * inv(partial Eq/partial x) * (partial Eq/partial bhp_target) + // where Eq(x)=0 is the well equation setup with bhp control and primary variables x // We shouldn't have zero rates at this stage, but check bool zero_rates; @@ -1363,6 +1363,7 @@ namespace Opm const EvalWell comp_rate = this->primary_variables_.getQs(comp_idx); const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx); for (size_t pvIdx = 0; pvIdx < num_eq; ++pvIdx) { + // well primary variable derivatives in EvalWell start at position Indices::numEq ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq); } ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value(); @@ -1610,7 +1611,9 @@ namespace Opm this->regularize_ = false; const auto& summary_state = ebosSimulator.vanguard().summaryState(); - // Max status switch frequency should be 2 to avoid getting stuck in cycle + // Always take a few (more than one) iterations after a switch before allowing a new switch + // The optimal number here is subject to further investigation, but it has been observerved + // that unless this number is >1, we may get stuck in a cycle const int min_its_after_switch = 3; int its_since_last_switch = min_its_after_switch; int switch_count= 0; @@ -1754,7 +1757,7 @@ namespace Opm } else { this->wellStatus_ = well_status_orig; this->operability_status_ = operability_orig; - const std::string message = fmt::format(" Well {} did not converged in {} inner iterations (" + const std::string message = fmt::format(" Well {} did not converge in {} inner iterations (" "{} control/status switches).", this->name(), it, switch_count); deferred_logger.debug(message); } diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 9a5ae4b88..6abd34b16 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -242,7 +242,7 @@ namespace Opm const double alq_value, DeferredLogger& deferred_logger) const override; - void updateIPRImplicit(const Simulator& ebosSimulator, WellState& well_state, DeferredLogger& deferred_logger); + void updateIPRImplicit(const Simulator& ebosSimulator, WellState& well_state, DeferredLogger& deferred_logger) override; virtual void computeWellRatesWithBhp( const Simulator& ebosSimulator, diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 3ea0da065..c74e87603 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -841,8 +841,8 @@ namespace Opm { // 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 + // dr/dbhp = - (partial r/partial x) * inv(partial Eq/partial x) * (partial Eq/partial bhp_target) + // where Eq(x)=0 is the well equation setup with bhp control and primary variables x // We shouldn't have zero rates at this stage, but check bool zero_rates; @@ -853,7 +853,7 @@ namespace Opm } auto& ws = well_state.well(this->index_of_well_); if (zero_rates) { - const auto msg = fmt::format("updateIPRImplicit: Well {} has zero rate, IPRs might be probelmatic", this->name()); + const auto msg = fmt::format("updateIPRImplicit: Well {} has zero rate, IPRs might be problematic", this->name()); deferred_logger.debug(msg); /* // could revert to standard approach here: @@ -882,7 +882,7 @@ namespace Opm 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(); + const size_t nEq = this->primary_variables_.numWellEq(); BVectorWell rhs(1); rhs[0].resize(nEq); // rhs = 0 except -1 for control eq @@ -899,6 +899,7 @@ namespace Opm EvalWell comp_rate = this->primary_variables_.getQs(comp_idx); const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx); for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) { + // well primary variable derivatives in EvalWell start at position Indices::numEq ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq); } ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value(); @@ -2321,14 +2322,15 @@ namespace Opm const bool fixed_status /*false*/) { 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 + // Always take a few (more than one) iterations after a switch before allowing a new switch + // The optimal number here is subject to further investigation, but it has been observerved + // that unless this number is >1, we may get stuck in a cycle constexpr int min_its_after_switch = 4; int its_since_last_switch = min_its_after_switch; int switch_count= 0; @@ -2341,7 +2343,7 @@ namespace Opm bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN); allow_switching = allow_switching && (!fixed_control || !fixed_status); bool changed = false; - bool final_check = false; + bool final_check = false; // well needs to be set operable or else solving/updating of re-opened wells is skipped this->operability_status_.resetOperability(); this->operability_status_.solvable = true; @@ -2364,7 +2366,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_) { @@ -2383,7 +2385,7 @@ namespace Opm its_since_last_switch = min_its_after_switch; } else { break; - } + } } ++it; @@ -2411,7 +2413,7 @@ namespace Opm } else { this->wellStatus_ = well_status_orig; this->operability_status_ = operability_orig; - const std::string message = fmt::format(" Well {} did not converged in {} inner iterations (" + const std::string message = fmt::format(" Well {} did not converge in {} inner iterations (" "{} 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 e1406610a..c7c847263 100644 --- a/opm/simulators/wells/VFPHelpers.cpp +++ b/opm/simulators/wells/VFPHelpers.cpp @@ -543,7 +543,7 @@ intersectWithIPR(const VFPProdTable& table, { // Given fixed thp, wfr, gfr and alq, this function finds a stable (-flo, bhp)-intersection // between the ipr-line and bhp(flo) if such an intersection exists. For multiple stable - // intersections, the one corresonding the largest flo i returned. + // intersections, the one corresponding the largest flo is returned. // The adjust_bhp-function is used to adjust the vfp-table bhp-values to actual bhp-values due // vfp/well ref-depth differences and/or WVFPDP-related pressure adjustments. diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index 938a2dc11..2d99866f3 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -962,7 +962,7 @@ getFloIPR(const WellState& well_state, const double& liquid_b = pu.phase_used[BlackoilPhases::Liquid]? ipr_b[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0; const double& vapour_b = pu.phase_used[BlackoilPhases::Vapour]? ipr_b[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0; // The getFlo helper is indended to pick one or add two of the phase rates (depending on FLO-type), - // but we can equally use it to pick/add the corresonding ipr_a, ipr_b + // but we can equally use it to pick/add the corresponding ipr_a, ipr_b return std::make_pair(detail::getFlo(table, aqua_a, liquid_a, vapour_a), detail::getFlo(table, aqua_b, liquid_b, vapour_b)); } diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 1f4676c0d..947a33c89 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -748,7 +748,7 @@ void WellInterfaceGeneric:: prepareForPotentialCalculations(const SummaryState& summary_state, WellState& well_state, Well::InjectionControls& inj_controls, - Well::ProductionControls& prod_controls) + Well::ProductionControls& prod_controls) const { const bool has_thp = this->wellHasTHPConstraints(summary_state); auto& ws = well_state.well(this->index_of_well_); diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 966336989..ca38a1352 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -248,7 +248,7 @@ protected: void prepareForPotentialCalculations(const SummaryState& summary_state, WellState& well_state, Well::InjectionControls& inj_controls, - Well::ProductionControls& prod_controls); + Well::ProductionControls& prod_controls) const; // definition of the struct OperabilityStatus struct OperabilityStatus { diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 51da204ca..db8a85d29 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -270,7 +270,7 @@ namespace Opm { 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; } @@ -278,7 +278,6 @@ namespace Opm const double sgn = this->isInjector() ? 1.0 : -1.0; if (!this->wellIsStopped()){ if (wqTotal*sgn <= 0.0 && !fixed_status){ - //std::cout << "Stopping well:" << this->name() << std::endl; this->stopWell(); return true; } else { @@ -322,13 +321,10 @@ namespace Opm // if the well can operate, it must at least be able to produce at the lowest bhp of the bhp-curve (explicit fractions) 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 = (this->isInjector())? inj_limit : prod_limit; if (has_thp) { @@ -357,9 +353,7 @@ 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(); @@ -367,7 +361,12 @@ namespace Opm initPrimaryVariablesEvaluation(); if (this->isProducer()) { - gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger); + 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); + } } WellTestState welltest_state_temp; @@ -462,7 +461,7 @@ namespace Opm 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); @@ -512,32 +511,29 @@ namespace Opm const bool isThp = ws.production_cmode == Well::ProducerCMode::THP; // check stability of solution under thp-control - if (true) { - if (converged && !this->stopppedOrZeroRateTarget(summary_state, well_state) && isThp) { - auto rates = well_state.well(this->index_of_well_).surface_rates; - this->adaptRatesForVFP(rates); - this->updateIPRImplicit(ebos_simulator, well_state, deferred_logger); - bool is_stable = WellBhpThpCalculator(*this).isStableSolution(well_state, this->well_ecl_, rates, summary_state); - if (!is_stable) { - // solution converged to an unstable point! - this->operability_status_.use_vfpexplicit = true; - // msg = ... - auto bhp_stable = WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state); - // if we find an intersection with a sufficiently lower bhp, re-solve equations - const double reltol = 1e-3; - const double cur_bhp = ws.bhp; - if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){ - const auto msg = fmt::format("isStableSolution: Well {} converged to an unstable solution, re-solving", this->name()); - deferred_logger.debug(msg); - solveWellWithBhp(ebos_simulator, dt, bhp_stable.value(), well_state, deferred_logger); - // re-solve with hopefully good initial guess - ws.thp = this->getTHPConstraint(summary_state); - converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); - } - } - } + if (converged && !this->stopppedOrZeroRateTarget(summary_state, well_state) && isThp) { + auto rates = well_state.well(this->index_of_well_).surface_rates; + this->adaptRatesForVFP(rates); + this->updateIPRImplicit(ebos_simulator, well_state, deferred_logger); + bool is_stable = WellBhpThpCalculator(*this).isStableSolution(well_state, this->well_ecl_, rates, summary_state); + 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); + // if we find an intersection with a sufficiently lower bhp, re-solve equations + const double reltol = 1e-3; + const double cur_bhp = ws.bhp; + if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){ + const auto msg = fmt::format("Well {} converged to an unstable solution, re-solving", this->name()); + deferred_logger.debug(msg); + solveWellWithBhp(ebos_simulator, dt, bhp_stable.value(), well_state, deferred_logger); + // re-solve with hopefully good initial guess + ws.thp = this->getTHPConstraint(summary_state); + converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + } + } } - + if (!converged) { // Well did not converge, switch to explicit fractions this->operability_status_.use_vfpexplicit = true; @@ -900,8 +896,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.has_well(well_name)) { - auto increment = glo.gaslift_increment(); + if(glo.active() && glo.has_well(well_name)) { + const auto increment = glo.gaslift_increment(); auto alq = well_state.getALQ(well_name); bool converged; while (alq > 0) { @@ -934,9 +930,8 @@ namespace Opm deferred_logger.info(msg); return; } - const auto& well_ecl = this->wellEcl(); const auto& schedule = ebos_simulator.vanguard().schedule(); - auto report_step_idx = ebos_simulator.episodeIndex(); + const 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( @@ -952,6 +947,7 @@ 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(); @@ -970,7 +966,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) { @@ -1016,7 +1012,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