OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4986 from steink/improved_local_thp_solve

Browse Source 
Implicit IPR for stability/operablity checks and solution of thp-wells
This commit is contained in:
Bård Skaflestad
2023-12-08 09:59:07 +01:00
committed by GitHub
parent 446670943a 051f4d4eab
commit 91fdded473
27 changed files with 1025 additions and 103 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
opm/simulators/flow/BlackoilModelParametersEbos.hpp
View File

@@ -182,6 +182,10 @@ template<class TypeTag, class MyTypeTag>
struct LocalWellSolveControlSwitching {
using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct UseImplicitIpr {
using type = UndefinedProperty;
};
// Network solver parameters
template<class TypeTag, class MyTypeTag>
struct NetworkMaxStrictIterations {
@@ -386,6 +390,10 @@ template<class TypeTag>
struct LocalWellSolveControlSwitching<TypeTag, TTag::FlowModelParameters> {
static constexpr bool value = false;
};
template<class TypeTag>
struct UseImplicitIpr<TypeTag, TTag::FlowModelParameters> {
static constexpr bool value = false;
};
// Network solver parameters
template<class TypeTag>
@@ -559,6 +567,9 @@ namespace Opm
/// Whether to allow control switching during local well solutions
bool local_well_solver_control_switching_;
/// Whether to use implicit IPR for thp stability checks and solution search
bool use_implicit_ipr_;
/// Maximum number of iterations in the network solver before relaxing tolerance
int network_max_strict_iterations_;
@@ -619,6 +630,7 @@ namespace Opm
max_number_of_well_switches_ = EWOMS_GET_PARAM(TypeTag, int, MaximumNumberOfWellSwitches);
use_average_density_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseAverageDensityMsWells);
local_well_solver_control_switching_ = EWOMS_GET_PARAM(TypeTag, bool, LocalWellSolveControlSwitching);
use_implicit_ipr_ = EWOMS_GET_PARAM(TypeTag, bool, UseImplicitIpr);
nonlinear_solver_ = EWOMS_GET_PARAM(TypeTag, std::string, NonlinearSolver);
const auto approach = EWOMS_GET_PARAM(TypeTag, std::string, LocalSolveApproach);
if (approach == "jacobi") {
@@ -680,6 +692,7 @@ namespace Opm
EWOMS_REGISTER_PARAM(TypeTag, int, MaximumNumberOfWellSwitches, "Maximum number of times a well can switch to the same control");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseAverageDensityMsWells, "Approximate segment densitities by averaging over segment and its outlet");
EWOMS_REGISTER_PARAM(TypeTag, bool, LocalWellSolveControlSwitching, "Allow control switching during local well solutions");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseImplicitIpr, "Compute implict IPR for stability checks and stable solution search");
EWOMS_REGISTER_PARAM(TypeTag, int, NetworkMaxStrictIterations, "Maximum iterations in network solver before relaxing tolerance");
EWOMS_REGISTER_PARAM(TypeTag, int, NetworkMaxIterations, "Maximum number of iterations in the network solver before giving up");
EWOMS_REGISTER_PARAM(TypeTag, std::string, NonlinearSolver, "Choose nonlinear solver. Valid choices are newton or nldd.");

10
opm/simulators/wells/MultisegmentWell.hpp
View File

@@ -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<double>& well_potentials,
DeferredLogger& deferred_logger) const;
virtual double getRefDensity() const override;
virtual bool iterateWellEqWithControl(const Simulator& ebosSimulator,
@@ -262,7 +268,9 @@ namespace Opm
const Well::ProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger) override;
DeferredLogger& deferred_logger,
const bool fixed_control = false,
const bool fixed_status = false) override;
virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
const double dt,

7
opm/simulators/wells/MultisegmentWellEquations.cpp
View File

@@ -181,6 +181,13 @@ MultisegmentWellEquations<Scalar,numWellEq,numEq>::solve() const
return mswellhelpers::applyUMFPack(*duneDSolver_, resWell_);
}
template<class Scalar, int numWellEq, int numEq>
typename MultisegmentWellEquations<Scalar,numWellEq,numEq>::BVectorWell
MultisegmentWellEquations<Scalar,numWellEq,numEq>::solve(const BVectorWell& rhs) const
{
return mswellhelpers::applyUMFPack(*duneDSolver_, rhs);
}
template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
recoverSolutionWell(const BVector& x, BVectorWell& xw) const

3
opm/simulators/wells/MultisegmentWellEquations.hpp
View File

@@ -96,6 +96,9 @@ public:
//! \brief Apply inverted D matrix to residual and return result.
BVectorWell solve() const;
//! \brief Apply inverted D matrix to rhs and return result.
BVectorWell solve(const BVectorWell& rhs) const;
//! \brief Recover well solution.
//! \details xw = inv(D)*(rw - C*x)
void recoverSolutionWell(const BVector& x, BVectorWell& xw) const;

9
opm/simulators/wells/MultisegmentWellEval.cpp
View File

@@ -86,7 +86,8 @@ getWellConvergence(const WellState& well_state,
const double relaxed_inner_tolerance_flow_ms_well,
const double tolerance_pressure_ms_wells,
const double relaxed_inner_tolerance_pressure_ms_well,
const bool relax_tolerance) const
const bool relax_tolerance,
const bool well_is_stopped) const
{
assert(int(B_avg.size()) == baseif_.numComponents());
@@ -161,12 +162,13 @@ getWellConvergence(const WellState& well_state,
tolerance_wells,
max_residual_allowed},
std::abs(linSys_.residual()[0][SPres]),
well_is_stopped,
report,
deferred_logger);
// for stopped well, we do not enforce the following checking to avoid dealing with sign of near-zero values
// for BHP or THP controlled wells, we need to make sure the flow direction is correct
if (!baseif_.wellIsStopped() && baseif_.isPressureControlled(well_state)) {
if (!well_is_stopped && baseif_.isPressureControlled(well_state)) {
// checking the flow direction
const double sign = baseif_.isProducer() ? -1. : 1.;
const auto weight_total_flux = this->primary_variables_.getWQTotal() * sign;
@@ -540,9 +542,6 @@ getResidualMeasureValue(const WellState& well_state,
++count;
}
// if (count == 0), it should be converged.
assert(count != 0);
return sum;
}

3
opm/simulators/wells/MultisegmentWellEval.hpp
View File

@@ -106,7 +106,8 @@ protected:
const double relaxed_inner_tolerance_flow_ms_well,
const double tolerance_pressure_ms_wells,
const double relaxed_inner_tolerance_pressure_ms_well,
const bool relax_tolerance) const;
const bool relax_tolerance,
const bool well_is_stopped) const;
std::pair<bool, std::vector<Scalar> >
getFiniteWellResiduals(const std::vector<Scalar>& B_avg,

203
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@@ -206,7 +206,8 @@ namespace Opm
this->param_.relaxed_tolerance_flow_well_,
this->param_.tolerance_pressure_ms_wells_,
this->param_.relaxed_tolerance_pressure_ms_well_,
relax_tolerance);
relax_tolerance,
this->wellIsStopped());
}
@@ -286,13 +287,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 +499,74 @@ namespace Opm
return potentials;
}
template<typename TypeTag>
bool
MultisegmentWell<TypeTag>::
computeWellPotentialsImplicit(const Simulator& ebos_simulator,
std::vector<double>& 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<TypeTag> 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
bool converged = false;
if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
converged = well_copy.solveWellWithTHPConstraint(ebos_simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
} else {
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 <typename TypeTag>
void
@@ -1227,6 +1305,73 @@ namespace Opm
}
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
updateIPRImplicit(const Simulator& ebos_simulator, 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 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;
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, IPRs might be problematic", this->name());
deferred_logger.debug(msg);
/*
// could revert to standard approach here:
updateIPR(ebos_simulator, deferred_logger);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
ws.implicit_ipr_a[idx] = this->ipr_a_[comp_idx];
ws.implicit_ipr_b[idx] = this->ipr_b_[comp_idx];
}
return;
*/
}
const auto& group_state = ebos_simulator.problem().wellModel().groupState();
std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
std::fill(ws.implicit_ipr_b.begin(), ws.implicit_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 = ebos_simulator.timeStepSize();
assembleWellEqWithoutIteration(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
BVectorWell rhs(this->numberOfSegments());
rhs = 0.0;
rhs[0][SPres] = -1.0;
const BVectorWell x_well = this->linSys_.solve(rhs);
constexpr int num_eq = MSWEval::numWellEq;
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
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();
}
// reset cmode
ws.production_cmode = cmode;
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
@@ -1441,10 +1586,10 @@ namespace Opm
const Well::ProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger)
DeferredLogger& deferred_logger,
const bool fixed_control /*false*/,
const bool fixed_status /*false*/)
{
//if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true;
const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
{
@@ -1464,25 +1609,40 @@ namespace Opm
[[maybe_unused]] int stagnate_count = 0;
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
const int min_its_after_switch = 2;
// 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;
const auto well_status = this->wellStatus_;
const auto& summary_state = ebosSimulator.vanguard().summaryState();
const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN);
// if we fail to solve eqs, we reset status/operability before leaving
const auto well_status_orig = this->wellStatus_;
const auto operability_orig = this->operability_status_;
auto well_status_cur = well_status_orig;
int status_switch_count = 0;
// only allow switcing if well is not under zero-rate target and is open from schedule
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;
// 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;
for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
its_since_last_switch++;
if (its_since_last_switch >= min_its_after_switch){
if (allow_switching && its_since_last_switch >= min_its_after_switch){
const double wqTotal = this->primary_variables_.getWQTotal().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, fixed_control, fixed_status);
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;
@@ -1586,9 +1746,6 @@ 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
// on the outside based on operability status
this->wellStatus_ = well_status;
}
std::string message = fmt::format(" Well {} converged in {} inner iterations ("
"{} control/status switches).", this->name(), it, switch_count);
@@ -1598,7 +1755,9 @@ namespace Opm
}
deferred_logger.debug(message);
} else {
const std::string message = fmt::format(" Well {} did not converged in {} inner iterations ("
this->wellStatus_ = well_status_orig;
this->operability_status_ = operability_orig;
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);
}

6
opm/simulators/wells/SingleWellState.cpp
View File

@@ -42,6 +42,8 @@ SingleWellState::SingleWellState(const std::string& name_,
, temperature(temp)
, well_potentials(pu_.num_phases)
, productivity_index(pu_.num_phases)
, implicit_ipr_a(pu_.num_phases)
, implicit_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->implicit_ipr_a.begin(), this->implicit_ipr_a.end(), 0);
std::fill(this->implicit_ipr_b.begin(), this->implicit_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->implicit_ipr_a == rhs.implicit_ipr_a &&
this->implicit_ipr_b == rhs.implicit_ipr_b &&
this->surface_rates == rhs.surface_rates &&
this->reservoir_rates == rhs.reservoir_rates &&
this->prev_surface_rates == rhs.prev_surface_rates &&

4
opm/simulators/wells/SingleWellState.hpp
View File

@@ -61,6 +61,8 @@ public:
serializer(phase_mixing_rates);
serializer(well_potentials);
serializer(productivity_index);
serializer(implicit_ipr_a);
serializer(implicit_ipr_b);
serializer(surface_rates);
serializer(reservoir_rates);
serializer(prev_surface_rates);
@@ -98,6 +100,8 @@ public:
std::vector<double> well_potentials;
std::vector<double> productivity_index;
std::vector<double> implicit_ipr_a;
std::vector<double> implicit_ipr_b;
std::vector<double> surface_rates;
std::vector<double> reservoir_rates;
std::vector<double> prev_surface_rates;

10
opm/simulators/wells/StandardWell.hpp
View File

@@ -212,7 +212,9 @@ namespace Opm
const Well::ProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger) override;
DeferredLogger& deferred_logger,
const bool fixed_control = false,
const bool fixed_status = false) override;
/// \brief Wether the Jacobian will also have well contributions in it.
virtual bool jacobianContainsWellContributions() const override
@@ -240,6 +242,8 @@ namespace Opm
const double alq_value,
DeferredLogger& deferred_logger) const override;
void updateIPRImplicit(const Simulator& ebosSimulator, WellState& well_state, DeferredLogger& deferred_logger) override;
virtual void computeWellRatesWithBhp(
const Simulator& ebosSimulator,
const double& bhp,
@@ -321,6 +325,10 @@ namespace Opm
DeferredLogger& deferred_logger,
const WellState &well_state) const;
bool computeWellPotentialsImplicit(const Simulator& ebos_simulator,
std::vector<double>& well_potentials,
DeferredLogger& deferred_logger) const;
virtual double getRefDensity() const override;
// get the mobility for specific perforation

5
opm/simulators/wells/StandardWellEquations.cpp
View File

@@ -177,6 +177,11 @@ void StandardWellEquations<Scalar,numEq>::solve(BVectorWell& dx_well) const
invDuneD_.mv(resWell_, dx_well);
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::solve(const BVectorWell& rhs_well, BVectorWell& x_well) const
{
invDuneD_.mv(rhs_well, x_well);
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::

3
opm/simulators/wells/StandardWellEquations.hpp
View File

@@ -89,6 +89,9 @@ public:
//! \brief Apply inverted D matrix to residual and store in vector.
void solve(BVectorWell& dx_well) const;
//! \brief Apply inverted D matrix to rhs and store in vector.
void solve(const BVectorWell& rhs_well, BVectorWell& x_well) const;
//! \brief Invert D matrix.
void invert();

4
opm/simulators/wells/StandardWellEval.cpp
View File

@@ -106,6 +106,7 @@ getWellConvergence(const WellState& well_state,
const double tol_wells,
const double relaxed_tolerance_flow,
const bool relax_tolerance,
const bool well_is_stopped,
std::vector<double>& res,
DeferredLogger& deferred_logger) const
{
@@ -150,12 +151,13 @@ getWellConvergence(const WellState& well_state,
checkConvergenceControlEq(well_state,
{1.e3, 1.e4, 1.e-4, 1.e-6, maxResidualAllowed},
std::abs(this->linSys_.residual()[0][Bhp]),
well_is_stopped,
report,
deferred_logger);
// for stopped well, we do not enforce the following checking to avoid dealing with sign of near-zero values
// for BHP or THP controlled wells, we need to make sure the flow direction is correct
if (!baseif_.wellIsStopped() && baseif_.isPressureControlled(well_state)) {
if (!well_is_stopped && baseif_.isPressureControlled(well_state)) {
// checking the flow direction
const double sign = baseif_.isProducer() ? -1. : 1.;
const auto weight_total_flux = this->primary_variables_.value(PrimaryVariables::WQTotal) * sign;

1
opm/simulators/wells/StandardWellEval.hpp
View File

@@ -85,6 +85,7 @@ protected:
const double tol_wells,
const double relaxed_tolerance_flow,
const bool relax_tolerance,
const bool well_is_stopped,
std::vector<double>& res,
DeferredLogger& deferred_logger) const;

232
opm/simulators/wells/StandardWell_impl.hpp
View File

@@ -832,6 +832,81 @@ namespace Opm
this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
}
template<typename TypeTag>
void
StandardWell<TypeTag>::
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 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;
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, IPRs might be problematic", this->name());
deferred_logger.debug(msg);
/*
// could revert to standard approach here:
updateIPR(ebos_simulator, deferred_logger);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
ws.implicit_ipr_a[idx] = this->ipr_a_[comp_idx];
ws.implicit_ipr_b[idx] = this->ipr_b_[comp_idx];
}
return;
*/
}
const auto& group_state = ebosSimulator.problem().wellModel().groupState();
std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
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 size_t 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) {
// 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();
}
// reset cmode
ws.production_cmode = cmode;
}
template<typename TypeTag>
void
@@ -1107,6 +1182,7 @@ namespace Opm
tol_wells,
this->param_.relaxed_tolerance_flow_well_,
relax_tolerance,
this->wellIsStopped(),
res,
deferred_logger);
@@ -1492,6 +1568,67 @@ namespace Opm
return potentials;
}
template<typename TypeTag>
bool
StandardWell<TypeTag>::
computeWellPotentialsImplicit(const Simulator& ebos_simulator,
std::vector<double>& 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<TypeTag> 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.
bool converged = false;
if (this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state)) {
converged = well_copy.solveWellWithTHPConstraint(ebos_simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger);
} else {
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<typename TypeTag>
@@ -1554,29 +1691,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<double>::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<double>::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 +2317,48 @@ namespace Opm
const Well::ProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger)
DeferredLogger& deferred_logger,
const bool fixed_control /*false*/,
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
constexpr int min_its_after_switch = 2;
// 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;
const auto well_status = this->wellStatus_;
const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN);
// if we fail to solve eqs, we reset status/operability before leaving
const auto well_status_orig = this->wellStatus_;
const auto operability_orig = this->operability_status_;
auto well_status_cur = well_status_orig;
int status_switch_count = 0;
// only allow switcing if well is not under zero-rate target and is open from schedule
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;
// 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;
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, fixed_control, fixed_status);
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;
@@ -2232,8 +2391,9 @@ namespace Opm
++it;
solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
initPrimaryVariablesEvaluation();
} while (it < max_iter);
} while (it < max_iter);
if (converged) {
if (allow_switching){
// update operability if status change
@@ -2244,18 +2404,12 @@ namespace Opm
} else {
this->operability_status_.operable_under_only_bhp_limit = !is_stopped;
}
// We reset the well status to its original state. Status is updated
// on the outside based on operability status
// \Note for future reference: For the well to update its status to stop/shut,
// the flag changed_to_stopped_this_step_ in prepareWellBeforeAssembling needs to be set to true.
// 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;
}
} else {
const std::string message = fmt::format(" Well {} did not converged in {} inner iterations ("
"{} control/status switches).", this->name(), it, switch_count);
this->wellStatus_ = well_status_orig;
this->operability_status_ = operability_orig;
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 ?
}

94
opm/simulators/wells/VFPHelpers.cpp
View File

@@ -499,6 +499,100 @@ double findTHP(const std::vector<double>& bhp_array,
return thp;
}
std::pair<double, double>
getMinimumBHPCoordinate(const VFPProdTable& table,
const double thp,
const double wfr,
const double gfr,
const double alq)
{
// Given fixed thp, wfr, gfr and alq, this function finds the minimum bhp and returns
// the corresponding pair (-flo_at_bhp_min, bhp_min). No assumption is taken on the
// shape of the function bhp(flo), so all points in the flo-axis is checked.
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;
const std::vector<double>& flos = table.getFloAxis();
for (size_t i = 0; i < flos.size(); ++i) {
flo_i = detail::findInterpData(flos[i], flos);
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<std::pair<double, double>>
intersectWithIPR(const VFPProdTable& table,
const double thp,
const double wfr,
const double gfr,
const double alq,
const double ipr_a,
const double ipr_b,
const std::function<double(const double)>& adjust_bhp)
{
// 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 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.
// 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_a, adjust_bhp(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 = adjust_bhp(bhp_i.value) - ipr_a/ipr_b; // +0.0/ipr_b
std::vector<double> 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 = adjust_bhp(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~0)
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 <typename T>
T getFlo(const VFPProdTable& table,
const T& aqua,

24
opm/simulators/wells/VFPHelpers.hpp
View File

@@ -25,6 +25,7 @@
#include <functional>
#include <map>
#include <vector>
#include <optional>
/**
* This file contains a set of helper functions used by VFPProd / VFPInj.
@@ -181,6 +182,29 @@ double findTHP(const std::vector<double>& bhp_array,
const std::vector<double>& thp_array,
double bhp);
/**
* Get (flo, bhp) at minimum bhp for given thp,wfr,gfr,alq
*/
std::pair<double, double>
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<std::pair<double, double>>
intersectWithIPR(const VFPProdTable& table,
const double thp,
const double wfr,
const double gfr,
const double alq,
const double ipr_a,
const double ipr_b,
const std::function<double(const double)>& adjust_bhp);
} // namespace detail

17
opm/simulators/wells/VFPProdProperties.cpp
View File

@@ -137,6 +137,22 @@ bhpwithflo(const std::vector<double>& 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 );
@@ -177,7 +193,6 @@ 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.setValue(bhp_val.value);
return bhp;
}

6
opm/simulators/wells/VFPProdProperties.hpp
View File

@@ -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<double> bhpwithflo(const std::vector<double>& flos,

128
opm/simulators/wells/WellBhpThpCalculator.cpp
View File

@@ -365,6 +365,7 @@ calculateBhpFromThp(const WellState& well_state,
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 bool use_vfpexplicit = well_.useVfpExplicit();
bhp_tab = well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number,
aqua, liquid, vapour,
thp_limit,
@@ -387,6 +388,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 +866,107 @@ bruteForceBracket(const std::function<double(const double)>& eq,
return bracket_found;
}
bool WellBhpThpCalculator::
isStableSolution(const WellState& well_state,
const Well& well,
const std::vector<double>& rates,
const SummaryState& summaryState) 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();
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 = getFloIPR(well_state, well, summaryState);
return bhp.dflo + 1/ipr.second >= 0;
}
}
std::optional<double> WellBhpThpCalculator::
estimateStableBhp(const WellState& well_state,
const Well& well,
const std::vector<double>& rates,
const double rho,
const SummaryState& summaryState) 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);
}
auto ipr = getFloIPR(well_state, well, summaryState);
const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth,
rho, well_.gravity());
auto bhp_adjusted = [this, &thp, &dp_hydro](const double bhp) {
return bhp - dp_hydro + getVfpBhpAdjustment(bhp, thp);
};
const auto retval = detail::intersectWithIPR(table, thp, wfr, gfr, well_.getALQ(well_state), ipr.first, ipr.second, bhp_adjusted);
if (retval.has_value()) {
// returned pair is (flo, bhp)
return retval.value().second;
} else {
return std::nullopt;
}
}
std::pair<double, double> WellBhpThpCalculator::
getFloIPR(const WellState& well_state,
const Well& well,
const SummaryState& summary_state) const
{
// Convert ipr_a's and ipr_b's to our particular choice of FLO
const auto& controls = well.productionControls(summary_state);
const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const auto& pu = well_.phaseUsage();
const auto& ipr_a = well_state.well(well_.indexOfWell()).implicit_ipr_a;
const double& aqua_a = pu.phase_used[BlackoilPhases::Aqua]? ipr_a[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0;
const double& liquid_a = pu.phase_used[BlackoilPhases::Liquid]? ipr_a[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0;
const double& vapour_a = pu.phase_used[BlackoilPhases::Vapour]? ipr_a[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0;
const auto& ipr_b = well_state.well(well_.indexOfWell()).implicit_ipr_b;
const double& aqua_b = pu.phase_used[BlackoilPhases::Aqua]? ipr_b[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0;
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 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));
}
#define INSTANCE(...) \
template __VA_ARGS__ WellBhpThpCalculator:: \
calculateBhpFromThp<__VA_ARGS__>(const WellState&, \

22
opm/simulators/wells/WellBhpThpCalculator.hpp
View File

@@ -98,6 +98,28 @@ 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<double>& rates,
const SummaryState& summaryState) const;
std::optional<double>
estimateStableBhp (const WellState& well_state,
const Well& well,
const std::vector<double>& rates,
const double rho,
const SummaryState& summaryState) const;
std::pair<double, double>
getFloIPR(const WellState& well_state,
const Well& well,
const SummaryState& summary_state) const;
private:
//! \brief Compute BHP from THP limit for an injector - implementation.
template<class ErrorPolicy>

3
opm/simulators/wells/WellConvergence.cpp
View File

@@ -37,6 +37,7 @@ void WellConvergence::
checkConvergenceControlEq(const WellState& well_state,
const Tolerances& tolerances,
const double well_control_residual,
const bool well_is_stopped,
ConvergenceReport& report,
DeferredLogger& deferred_logger) const
{
@@ -46,7 +47,7 @@ checkConvergenceControlEq(const WellState& well_state,
const int well_index = well_.indexOfWell();
const auto& ws = well_state.well(well_index);
if (well_.wellIsStopped()) {
if (well_is_stopped) {
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = tolerances.rates; // use smaller tolerance for zero control?
}

1
opm/simulators/wells/WellConvergence.hpp
View File

@@ -52,6 +52,7 @@ public:
void checkConvergenceControlEq(const WellState& well_state,
const Tolerances& tolerances,
const double well_control_residual,
const bool well_is_stopped,
ConvergenceReport& report,
DeferredLogger& deferred_logger) const;

37
opm/simulators/wells/WellInterface.hpp
View File

@@ -247,7 +247,9 @@ public:
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
const double WQTotal,
DeferredLogger& deferred_logger);
DeferredLogger& deferred_logger,
const bool fixed_control = false,
const bool fixed_status = false);
virtual void updatePrimaryVariables(const SummaryState& summary_state,
const WellState& well_state,
@@ -414,7 +416,13 @@ protected:
const WellProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger) = 0;
DeferredLogger& deferred_logger,
const bool fixed_control = false,
const bool fixed_status = false) = 0;
virtual void updateIPRImplicit(const Simulator& ebosSimulator,
WellState& well_state,
DeferredLogger& deferred_logger) = 0;
bool iterateWellEquations(const Simulator& ebosSimulator,
const double dt,
@@ -422,6 +430,31 @@ protected:
const GroupState& group_state,
DeferredLogger& deferred_logger);
bool solveWellWithTHPConstraint(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<double> estimateOperableBhp(const Simulator& ebos_simulator,
const double dt,
WellState& well_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 solveWellWithZeroRate(const Simulator& ebos_simulator,
const double dt,
WellState& well_state,
DeferredLogger& deferred_logger);
bool solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const GroupState& group_state,
DeferredLogger& deferred_logger);

43
opm/simulators/wells/WellInterfaceGeneric.cpp
View File

@@ -744,4 +744,47 @@ checkNegativeWellPotentials(std::vector<double>& well_potentials,
}
}
void WellInterfaceGeneric::
prepareForPotentialCalculations(const SummaryState& summary_state,
WellState& well_state,
Well::InjectionControls& inj_controls,
Well::ProductionControls& prod_controls) const
{
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

5
opm/simulators/wells/WellInterfaceGeneric.hpp
View File

@@ -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) const;
// definition of the struct OperabilityStatus
struct OperabilityStatus {
bool isOperableAndSolvable() const {

235
opm/simulators/wells/WellInterface_impl.hpp
View File

@@ -264,7 +264,9 @@ namespace Opm
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
const double wqTotal,
DeferredLogger& deferred_logger)
DeferredLogger& deferred_logger,
const bool fixed_control,
const bool fixed_status)
{
const auto& summary_state = ebos_simulator.vanguard().summaryState();
const auto& schedule = ebos_simulator.vanguard().schedule();
@@ -275,60 +277,65 @@ namespace Opm
const double sgn = this->isInjector() ? 1.0 : -1.0;
if (!this->wellIsStopped()){
if (wqTotal*sgn <= 0.0){
if (wqTotal*sgn <= 0.0 && !fixed_status){
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 (!fixed_control) {
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);
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;
}
} else {
} else if (!fixed_status){
// well is stopped, check if current bhp allows reopening
const double bhp = well_state.well(this->index_of_well_).bhp;
double prod_limit = prod_controls.bhp_limit;
double inj_limit = inj_controls.bhp_limit;
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.
std::vector<double> rates(this->num_components_);
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);
// 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);
}
}
const double bhp_diff = (this->isInjector())? inj_limit - bhp: bhp - prod_limit;
if (bhp_diff > 0){
this->openWell();
well_state.well(this->index_of_well_).bhp = (this->isInjector())? inj_limit : prod_limit;
if (has_thp) {
well_state.well(this->index_of_well_).thp = this->getTHPConstraint(summary_state);
}
return true;
} else {
return false;
}
} else {
return false;
}
}
@@ -448,7 +455,11 @@ 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->param_.use_implicit_ipr_ && this->well_ecl_.isProducer() && this->wellHasTHPConstraints(summary_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN)) {
converged = solveWellWithTHPConstraint(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 ) {
@@ -459,6 +470,174 @@ namespace Opm
return converged;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
solveWellWithTHPConstraint(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 = true;
auto& ws = well_state.well(this->index_of_well_);
// if well is stopped, check if we can reopen
if (this->wellIsStopped()) {
this->openWell();
auto bhp_target = estimateOperableBhp(ebos_simulator, dt, well_state, summary_state, deferred_logger);
if (!bhp_target.has_value()) {
// no intersection with ipr
const auto msg = fmt::format("estimateOperableBhp: Did not find operable BHP for well {}", this->name());
deferred_logger.debug(msg);
is_operable = false;
// solve with zero rates
solveWellWithZeroRate(ebos_simulator, dt, well_state, deferred_logger);
this->stopWell();
} else {
// solve well with the estimated target bhp (or limit)
ws.thp = this->getTHPConstraint(summary_state);
const double bhp = std::max(bhp_target.value(), prod_controls.bhp_limit);
solveWellWithBhp(ebos_simulator, dt, bhp, well_state, deferred_logger);
}
}
// solve well-equation
if (is_operable) {
converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
}
const bool isThp = ws.production_cmode == Well::ProducerCMode::THP;
// check stability of solution under thp-control
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;
this->openWell();
auto bhp_target = estimateOperableBhp(ebos_simulator, dt, well_state, summary_state, deferred_logger);
if (!bhp_target.has_value()) {
// well can't operate using explicit fractions
is_operable = false;
// solve with zero rate
converged = solveWellWithZeroRate(ebos_simulator, dt, well_state, deferred_logger);
this->stopWell();
} else {
// solve well with the estimated target bhp (or limit)
const double bhp = std::max(bhp_target.value(), prod_controls.bhp_limit);
solveWellWithBhp(ebos_simulator, dt, bhp, well_state, deferred_logger);
ws.thp = this->getTHPConstraint(summary_state);
converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
}
}
// update operability
is_operable = is_operable && !this->wellIsStopped();
this->operability_status_.can_obtain_bhp_with_thp_limit = is_operable;
this->operability_status_.obey_thp_limit_under_bhp_limit = is_operable;
return converged;
}
template<typename TypeTag>
std::optional<double>
WellInterface<TypeTag>::
estimateOperableBhp(const Simulator& ebos_simulator,
const double dt,
WellState& well_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 || 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);
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
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, /*fixed_control*/true);
ws.injection_cmode = cmode_inj;
ws.production_cmode = cmode_prod;
return converged;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
solveWellWithZeroRate(const Simulator& ebos_simulator,
const double dt,
WellState& well_state,
DeferredLogger& deferred_logger)
{
// Solve a well as stopped
const auto well_status_orig = this->wellStatus_;
this->stopWell();
auto group_state = GroupState(); // empty group
auto inj_controls = Well::InjectionControls(0);
auto prod_controls = Well::ProductionControls(0);
const bool converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger, /*fixed_control*/true, /*fixed_status*/ true);
this->wellStatus_ = well_status_orig;
return converged;
}
template<typename TypeTag>
bool