This commit is contained in:
Halvor Møll Nilsen 2025-02-13 14:49:34 +01:00 committed by GitHub
commit d7ca7a473e
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
21 changed files with 187 additions and 18 deletions

View File

@ -1240,6 +1240,7 @@ namespace Opm {
OPM_BEGIN_PARALLEL_TRY_CATCH();
{
// Set the well primary variables based on the value of well solutions
updatePrimaryVariables(local_deferredLogger);
initPrimaryVariablesEvaluation();
alq_updated = gaslift_.maybeDoGasLiftOptimize(simulator_,

View File

@ -28,6 +28,7 @@
#include <opm/simulators/wells/WellGroupHelpers.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <opm/common/TimingMacros.hpp>
#include <cassert>
namespace Opm::WGHelpers {
@ -116,6 +117,7 @@ guideRateSum(const Group& group,
const std::string& always_included_child,
const bool always_use_potentials)
{
OPM_TIMEFUNCTION();
Scalar total_guide_rate = 0.0;
int number_of_included_well_or_groups = 0;
for (const std::string& child_group : group.groups()) {
@ -159,6 +161,7 @@ guideRate(const std::string& name,
const std::string& always_included_child,
const bool always_use_potentials)
{
OPM_TIMEFUNCTION();
if (schedule_.hasWell(name, report_step_)) {
return WellGroupHelpers<Scalar>::getGuideRate(name, schedule_, well_state_, group_state_,
report_step_, guide_rate_, target_, pu_);

View File

@ -55,6 +55,7 @@ public:
const WellInterfaceGeneric<Scalar>& getWell() const override { return well_; }
private:
void solveWellWithTHPConstraintAlqImplicit(WellState<Scalar>& wellState,const GroupState<Scalar>& groupState) const override;
std::optional<Scalar>
computeBhpAtThpLimit_(Scalar alq,
bool debug_ouput = true) const override;

View File

@ -29,6 +29,7 @@
#include <opm/simulators/wells/GasLiftWellState.hpp>
#include <opm/simulators/wells/GroupState.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <fmt/format.h>
@ -420,6 +421,7 @@ std::optional<typename GasLiftSingleWellGeneric<Scalar>::LimitedRates>
GasLiftSingleWellGeneric<Scalar>::
computeLimitedWellRatesWithALQ_(Scalar alq) const
{
OPM_TIMEFUNCTION();
std::optional<LimitedRates> limited_rates;
if (auto rates = computeWellRatesWithALQ_(alq); rates) {
limited_rates = getLimitedRatesFromRates_(*rates);
@ -433,13 +435,41 @@ GasLiftSingleWellGeneric<Scalar>::
computeWellRatesWithALQ_(Scalar alq) const
{
OPM_TIMEFUNCTION();
std::optional<BasicRates> rates;
auto bhp_opt = computeBhpAtThpLimit_(alq);
if (bhp_opt) {
auto [bhp, bhp_is_limited] = getBhpWithLimit_(*bhp_opt);
rates = computeWellRates_(bhp, bhp_is_limited);
bool new_code = true;
if(new_code){
std::optional<BasicRates> rates_opt;
auto copy_well_state(this->well_state_);
copy_well_state.setALQ(this->well_name_, alq);
auto& ws = copy_well_state.well(this->well_name_);
//bool solved =
solveWellWithTHPConstraintAlqImplicit(copy_well_state,this->group_state_);
auto bhp_tmp = ws.bhp;
auto [bhp, bhp_is_limited] = getBhpWithLimit_(bhp_tmp);
bool tmp = bhp_is_limited;
BasicRates rates(-ws.surface_rates[this->oil_pos_],
-ws.surface_rates[this->gas_pos_],
-ws.surface_rates[this->water_pos_],
tmp);
if(ws.status == WellStatus::OPEN){
rates_opt = rates;
}
return rates_opt;
}else{
// old code
std::optional<BasicRates> rates;
auto bhp_opt = computeBhpAtThpLimit_(alq);
if (bhp_opt) {
auto [bhp, bhp_is_limited] = getBhpWithLimit_(*bhp_opt);
rates = computeWellRates_(bhp, bhp_is_limited);
}
return rates;
}
return rates;
//}
}
template<class Scalar>

View File

@ -152,7 +152,19 @@ protected:
bhp_is_limited = rates.bhp_is_limited;
return *this;
}
bool operator==(const BasicRates& rates){
return oil == rates.oil &&
gas == rates.gas &&
water == rates.water &&
bhp_is_limited == rates.bhp_is_limited;
}
bool approx(const BasicRates& rates, double tol){
return std::abs(oil - rates.oil) < tol &&
std::abs(gas - rates.gas) < tol &&
std::abs(water - rates.water) < tol &&
bhp_is_limited == rates.bhp_is_limited;
}
// This copy constructor cannot be defined inline here since LimitedRates
// has not been defined yet (it is defined below). Instead it is defined in
// in the .cpp file
@ -285,6 +297,7 @@ protected:
bool checkInitialALQmodified_(Scalar alq, Scalar initial_alq) const;
virtual bool checkThpControl_() const = 0;
virtual void solveWellWithTHPConstraintAlqImplicit(WellState<Scalar>& wellState,const GroupState<Scalar>& groupState) const = 0;
virtual std::optional<Scalar > computeBhpAtThpLimit_(Scalar alq,
bool debug_output = true) const = 0;

View File

@ -138,6 +138,19 @@ computeWellRates_(Scalar bhp, bool bhp_is_limited, bool debug_output ) const
};
}
template<typename TypeTag>
void GasLiftSingleWell<TypeTag>::solveWellWithTHPConstraintAlqImplicit(WellState<Scalar>& wellState,
const GroupState<Scalar>& groupState) const
{
OPM_TIMEFUNCTION();
this->well_.solveWellWithTHPConstraintConst(this->simulator_,wellState,groupState);
#if 0
auto well_copy(this->well_);
//const_cast<WellInterface<TypeTag>&>(this->well_).solveWellWithTHPConstraintALQ(this->simulator_,wellState,groupState);
this->well_.solveWellWithTHPConstraintALQ(this->simulator_,wellState,groupState);
#endif
}
template<typename TypeTag>
std::optional<typename GasLiftSingleWell<TypeTag>::Scalar>
GasLiftSingleWell<TypeTag>::

View File

@ -28,6 +28,7 @@
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/simulators/wells/GroupState.hpp>
#include <opm/common/TimingMacros.hpp>
namespace Opm {
@ -337,6 +338,7 @@ template<class Scalar>
bool GroupState<Scalar>::
has_production_control(const std::string& gname) const
{
OPM_TIMEFUNCTION();
auto group_iter = this->production_controls.find(gname);
if (group_iter == this->production_controls.end())
return false;
@ -356,6 +358,7 @@ template<class Scalar>
Group::ProductionCMode
GroupState<Scalar>::production_control(const std::string& gname) const
{
OPM_TIMEFUNCTION();
auto group_iter = this->production_controls.find(gname);
if (group_iter == this->production_controls.end())
throw std::logic_error("Could not find any control for production group: " + gname);

View File

@ -165,7 +165,11 @@ namespace Opm {
std::vector<Scalar> getPrimaryVars() const override;
int setPrimaryVars(typename std::vector<Scalar>::const_iterator it) override;
void solveWellWithTHPConstraintConst(const Simulator& simulator, WellState<Scalar>& wellState,
const GroupState<Scalar>& groupState) const{
auto copy_well(*this);
copy_well.solveWellWithTHPConstraintALQ(simulator, wellState, groupState);
}
protected:
// regularize msw equation
bool regularize_;

View File

@ -1475,6 +1475,7 @@ namespace Opm
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger)
{
OPM_TIMEFUNCTION();
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true;
const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
@ -1628,6 +1629,7 @@ namespace Opm
const bool fixed_control /*false*/,
const bool fixed_status /*false*/)
{
OPM_TIMEFUNCTION();
const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
{

View File

@ -245,13 +245,13 @@ update_producer_targets(const Well& ecl_well, const SummaryState& st)
case Well::ProducerCMode::THP:
case Well::ProducerCMode::BHP:
if (this->pu.phase_used[BlackoilPhases::Liquid]) {
this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = -1000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = 0.0 * (-1000.0) * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
}
if (this->pu.phase_used[BlackoilPhases::Aqua]) {
this->surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = -1000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
this->surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = 0.0 * (-1000.0) * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
}
if (this->pu.phase_used[BlackoilPhases::Vapour]){
this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = -100000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = 0.0 * (-100000.0) * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
}
break;

View File

@ -259,6 +259,21 @@ namespace Opm
int setPrimaryVars(typename std::vector<Scalar>::const_iterator it) override;
void solveWellWithTHPConstraintConst(const Simulator& simulator, WellState<Scalar>& wellState,
const GroupState<Scalar>& groupState) const{
OPM_TIMEFUNCTION();
auto well_copy(*this);
//well_copy.prepareForPotentialCalculations(summaryw_state, well_state_copy, inj_controls, prod_controls);
DeferredLogger deferred_logger;
//well_copy.calculateExplicitQuantities(simulator, wellState, deferred_logger);
well_copy.updatePrimaryVariables(simulator, wellState, deferred_logger);
well_copy.initPrimaryVariablesEvaluation();
well_copy.solveWellWithTHPConstraintALQ(simulator, wellState, groupState);
auto& ws = wellState.well(this->name());
if(ws.status == Opm::WellStatus::OPEN){
assert(!well_copy.wellIsStopped());
}
}
protected:
bool regularize_;

View File

@ -223,6 +223,7 @@ update(const WellState<Scalar>& well_state,
// BHP
value_[Bhp] = ws.bhp;
init();
}
template<class FluidSystem, class Indices>
@ -239,6 +240,7 @@ updatePolyMW(const WellState<Scalar>& well_state)
value_[Bhp + 1 + well_.numPerfs() + perf] = skin_pressure[perf];
}
}
init();
}
template<class FluidSystem, class Indices>
@ -301,6 +303,7 @@ updateNewton(const BVectorWell& dwells,
// so that bhp constaint can be an active control when needed.
constexpr Scalar bhp_lower_limit = 1. * unit::barsa - 1. * unit::Pascal;
value_[Bhp] = std::max(value_[Bhp] - dx1_limited, bhp_lower_limit);
init();
}
template<class FluidSystem, class Indices>
@ -320,6 +323,7 @@ updateNewtonPolyMW(const BVectorWell& dwells)
value_[pskin_index] -= relaxation_factor * dx_pskin;
}
}
init();
}
template<class FluidSystem, class Indices>

View File

@ -140,11 +140,11 @@ public:
//! \brief Returns a value.
Scalar value(const int idx) const
{ return value_[idx]; }
{ assert(evaluation_[idx].value() == value_[idx]);return value_[idx]; }
//! \brief Returns a const ref to an evaluation.
const EvalWell& eval(const int idx) const
{ return evaluation_[idx]; }
{ assert(evaluation_[idx].value() == value_[idx]); return evaluation_[idx]; }
//! \brief Set a value. Note that this does not also set the corresponding evaluation.
void setValue(const int idx, const Scalar val)

View File

@ -356,6 +356,7 @@ namespace Opm
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger)
{
OPM_TIMEFUNCTION();
// TODO: only_wells should be put back to save some computation
// for example, the matrices B C does not need to update if only_wells
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
@ -1751,6 +1752,7 @@ namespace Opm
std::vector<Scalar>& well_potentials,
DeferredLogger& deferred_logger) // const
{
OPM_TIMEFUNCTION();
const auto [compute_potential, bhp_controlled_well] =
this->WellInterfaceGeneric<Scalar>::computeWellPotentials(well_potentials, well_state);
@ -2354,6 +2356,7 @@ namespace Opm
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger)
{
OPM_TIMEFUNCTION();
const int max_iter = this->param_.max_inner_iter_wells_;
int it = 0;
bool converged;
@ -2402,6 +2405,7 @@ namespace Opm
const bool fixed_control /*false*/,
const bool fixed_status /*false*/)
{
OPM_TIMEFUNCTION();
const int max_iter = this->param_.max_inner_iter_wells_;
int it = 0;
bool converged = false;

View File

@ -29,6 +29,7 @@
#include <cassert>
#include <stdexcept>
#include <opm/common/TimingMacros.hpp>
namespace Opm::WGHelpers {
template<class Scalar>
@ -55,6 +56,7 @@ template<class Scalar>
template <typename RateType>
RateType TargetCalculator<Scalar>::calcModeRateFromRates(const RateType* rates) const
{
OPM_TIMEFUNCTION();
switch (cmode_) {
case Group::ProductionCMode::ORAT: {
assert(pu_.phase_used[BlackoilPhases::Liquid]);
@ -97,6 +99,7 @@ Scalar TargetCalculator<Scalar>::
groupTarget(const std::optional<Group::ProductionControls>& ctrl,
DeferredLogger& deferred_logger) const
{
OPM_TIMEFUNCTION();
if (!ctrl && !use_gpmaint_) {
OPM_DEFLOG_THROW(std::logic_error,
"Production group " + this->group_name_
@ -174,6 +177,7 @@ InjectionTargetCalculator(const Group::InjectionCMode& cmode,
, use_gpmaint_(use_gpmaint)
{
OPM_TIMEFUNCTION();
// initialize to avoid warning
pos_ = pu.phase_pos[BlackoilPhases::Aqua];
target_ = GuideRateModel::Target::WAT;
@ -206,6 +210,7 @@ Scalar InjectionTargetCalculator<Scalar>::
groupTarget(const std::optional<Group::InjectionControls>& ctrl,
DeferredLogger& deferred_logger) const
{
OPM_TIMEFUNCTION();
if (!ctrl && !use_gpmaint_) {
OPM_DEFLOG_THROW(std::logic_error,
"Injection group " + this->group_name_

View File

@ -644,6 +644,7 @@ bhpMax(const std::function<Scalar(const Scalar)>& fflo,
const Scalar vfp_flo_front,
DeferredLogger& deferred_logger) const
{
OPM_TIMEFUNCTION();
// Find the bhp-point where production becomes nonzero.
Scalar bhp_max = 0.0;
Scalar low = bhp_limit;
@ -724,6 +725,7 @@ computeBhpAtThpLimit(const std::function<std::vector<Scalar>(const Scalar)>& fra
const std::array<Scalar, 2>& range,
DeferredLogger& deferred_logger) const
{
//OPM_TIMEFUNCTION();
// Given a VFP function returning bhp as a function of phase
// rates and thp:
// fbhp(rates, thp),
@ -741,6 +743,7 @@ computeBhpAtThpLimit(const std::function<std::vector<Scalar>(const Scalar)>& fra
// highest rate) should be returned.
// Define the equation we want to solve.
auto eq = [&fbhp, &frates](Scalar bhp) {
return fbhp(frates(bhp)) - bhp;
};
@ -775,6 +778,7 @@ computeBhpAtThpLimit(const std::function<std::vector<Scalar>(const Scalar)>& fra
const Scalar bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0;
try {
OPM_TIMEBLOCK(RegularFalsiBisection);
const Scalar solved_bhp = RegulaFalsiBisection<ThrowOnError>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
return solved_bhp;
@ -794,6 +798,7 @@ bisectBracket(const std::function<Scalar(const Scalar)>& eq,
std::optional<Scalar>& approximate_solution,
DeferredLogger& deferred_logger) const
{
OPM_TIMEFUNCTION();
bool finding_bracket = false;
low = range[0];
high = range[1];
@ -869,6 +874,7 @@ bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
Scalar& low, Scalar& high,
DeferredLogger& deferred_logger)
{
OPM_TIMEFUNCTION();
bool bracket_found = false;
low = range[0];
high = range[1];
@ -901,6 +907,7 @@ isStableSolution(const WellState<Scalar>& well_state,
const std::vector<Scalar>& rates,
const SummaryState& summaryState) const
{
OPM_TIMEFUNCTION();
assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
assert(well_.isProducer()); // only valid for producers
@ -1013,6 +1020,7 @@ bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
const Scalar& limit,
DeferredLogger& deferred_logger)
{
OPM_TIMEFUNCTION();
bool bracket_found = false;
low = range[0];
high = range[1];
@ -1049,6 +1057,7 @@ WellBhpThpCalculator<Scalar>::
bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
Scalar& min_thp, Scalar& max_thp)
{
OPM_TIMEFUNCTION();
bool bracket_found = false;
constexpr int sample_number = 1000;
constexpr Scalar interval = 1E5;

View File

@ -27,6 +27,7 @@
#include <unordered_map>
#include <vector>
#include <opm/common/TimingMacros.hpp>
namespace Opm {
@ -151,6 +152,7 @@ public:
}
std::optional<int> well_index(const std::string& wname) const {
OPM_TIMEFUNCTION();
auto index_iter = this->index_map.find(wname);
if (index_iter != this->index_map.end())
return index_iter->second;
@ -159,6 +161,7 @@ public:
}
const std::string& well_name(std::size_t well_index) const {
OPM_TIMEFUNCTION();
for (const auto& [wname, windex] : this->index_map) {
if (windex == well_index)
return wname;

View File

@ -92,7 +92,7 @@ namespace Opm {
const int phasePos,
const bool injector)
{
OPM_TIMEFUNCTION();
Scalar rate = 0.0;
for (const std::string& groupName : group.groups()) {
const auto& groupTmp = schedule.getGroup(groupName, reportStepIdx);
@ -102,6 +102,7 @@ namespace Opm {
// only sum satelite production once
if (wellState.isRank0() && !injector) {
OPM_TIMEBLOCK(sumWellPhaseRatesSatelite);
const auto rateComp = selectRateComponent(wellState.phaseUsage(), phasePos);
if (rateComp.has_value()) {
rate += satelliteProduction(schedule[reportStepIdx], group.groups(), *rateComp);
@ -1062,6 +1063,7 @@ getGuideRate(const std::string& name,
const GuideRateModel::Target target,
const PhaseUsage& pu)
{
OPM_TIMEFUNCTION();
if (schedule.hasWell(name, reportStepIdx)) {
if (guideRate->has(name) || guideRate->hasPotentials(name)) {
return guideRate->get(name, target, getWellRateVector(wellState, pu, name));
@ -1174,6 +1176,8 @@ groupControlledWells(const Schedule& schedule,
OPM_TIMEFUNCTION();
const Group& group = schedule.getGroup(group_name, report_step);
int num_wells = 0;
{
OPM_TIMEBLOCK(iterateGroups);
for (const std::string& child_group : group.groups()) {
bool included = (child_group == always_included_child);
@ -1190,6 +1194,9 @@ groupControlledWells(const Schedule& schedule,
+= groupControlledWells(schedule, well_state, group_state, summary_state, guideRate, report_step, child_group, always_included_child, is_production_group, injection_phase);
}
}
}
{
OPM_TIMEBLOCK(iterateWells);
for (const std::string& child_well : group.wells()) {
bool included = (child_well == always_included_child);
if (is_production_group) {
@ -1199,6 +1206,7 @@ groupControlledWells(const Schedule& schedule,
}
const auto ctrl1 = group_state.production_control(group.name());
if (group.as_choke() && ((ctrl1 == Group::ProductionCMode::FLD) || (ctrl1 == Group::ProductionCMode::NONE))){
OPM_TIMEBLOCK(chokeCalculations);
// The auto choke group has not own group control but inherits control from an ancestor group.
// Number of wells should be calculated as zero when wells of auto choke group do not deliver target.
// This behaviour is then similar to no-autochoke group with wells not on GRUP control.
@ -1266,6 +1274,7 @@ groupControlledWells(const Schedule& schedule,
++num_wells;
}
}
}
return num_wells;
}
@ -1457,6 +1466,7 @@ checkGroupConstraintsProd(const std::string& name,
}
// check whether guide rate is violated
for (std::size_t ii = 1; ii < num_ancestors; ++ii) {
OPM_TIMEBLOCK(checkGuideRateVialted);
if (guideRate->has(chain[ii])) {
const auto& guided_group = chain[ii];
const Scalar grefficiency
@ -1473,6 +1483,7 @@ checkGroupConstraintsProd(const std::string& name,
// Compute portion of target corresponding to current_rate_available
Scalar target = orig_target;
for (std::size_t ii = 0; ii < num_ancestors; ++ii) {
OPM_TIMEBLOCK(computeTargetAvailable);
if ((ii == 0) || guideRate->has(chain[ii])) {
// Apply local reductions only at the control level
// (top) and for levels where we have a specified

View File

@ -339,6 +339,10 @@ public:
computeCurrentWellRates(const Simulator& simulator,
DeferredLogger& deferred_logger) const = 0;
virtual void solveWellWithTHPConstraintConst(const Simulator& simulator, WellState<Scalar>& wellState,
const GroupState<Scalar>& groupState) const = 0;
/// Modify the well_state's rates if there is only one nonzero rate.
/// If so, that rate is kept as is, but the others are set proportionally
/// to the rates returned by computeCurrentWellRates().
@ -439,6 +443,10 @@ protected:
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger);
bool solveWellWithTHPConstraintALQ(const Simulator& simulator, WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state);
bool solveWellWithTHPConstraint(const Simulator& simulator,
const double dt,
const Well::InjectionControls& inj_controls,
@ -506,7 +514,7 @@ protected:
const SingleWellState<Scalar>& ws) const;
};
} // namespace Opm
} // namespace OpmS
#include "WellInterface_impl.hpp"

View File

@ -295,7 +295,16 @@ namespace Opm
}
const bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
if (oscillating || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) || !(this->well_ecl_.getStatus() == WellStatus::OPEN)) {
// do not check groups zero rate in second case
bool zero_target = false;
if(this->param_.check_group_constraints_inner_well_iterations_){
zero_target = this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger);
}else{
const auto& summaryState = simulator.vanguard().summaryState();
zero_target = this->wellUnderZeroRateTargetIndividual(summaryState, well_state);
}
if (oscillating || zero_target || !(this->well_ecl_.getStatus() == WellStatus::OPEN)) {
return false;
}
@ -396,7 +405,7 @@ namespace Opm
WellState<Scalar> well_state_copy = well_state;
auto& ws = well_state_copy.well(this->indexOfWell());
calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger);
calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
updatePrimaryVariables(simulator, well_state_copy, deferred_logger);
@ -527,6 +536,30 @@ namespace Opm
return converged;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
solveWellWithTHPConstraintALQ(const Simulator& simulator, WellState<Scalar>& well_state, const GroupState<Scalar>& group_state){
OPM_TIMEFUNCTION();
const auto& summary_state = simulator.vanguard().summaryState();
auto inj_controls = this->well_ecl_.isInjector() ? this->well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
auto prod_controls = this->well_ecl_.isProducer() ? this->well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
const double dt = simulator.timeStepSize();
DeferredLogger deferred_logger;
this->prepareForPotentialCalculations(summary_state, well_state, inj_controls, prod_controls);
bool converged = this->solveWellWithTHPConstraint(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
auto& ws = well_state.well(this->name());
if(this->wellIsStopped()){
ws.status = WellStatus::STOP;
}else{
ws.status = WellStatus::OPEN;
}
assert(converged);
return converged;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
@ -845,6 +878,7 @@ namespace Opm
DeferredLogger& deferred_logger)
{
OPM_TIMEFUNCTION();
initPrimaryVariablesEvaluation();
const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
if (this->param_.check_well_operability_iter_)
@ -1547,6 +1581,7 @@ namespace Opm
DeferredLogger& deferred_logger,
const std::optional<bool> group_control) const
{
OPM_TIMEFUNCTION();
// Check if well is under zero rate target from group
const bool isGroupControlled = group_control.value_or(this->wellUnderGroupControl(well_state.well(this->index_of_well_)));
if (isGroupControlled) {
@ -1941,9 +1976,11 @@ namespace Opm
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const
{
OPM_TIMEFUNCTION();
OPM_TIMEFUNCTION();
const auto& summary_state = simulator.vanguard().summaryState();
// if(std::abs(this->getTHPConstraint(summary_state)-ws.thp) < 10*bar) {
// return true;
// }
auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
simulator, summary_state, this->getALQ(well_state), deferred_logger, /*iterate_if_no_solution */ true);
if (bhp_at_thp_limit) {

View File

@ -25,6 +25,7 @@
#include <dune/common/parallel/mpihelper.hh>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/TimingMacros.hpp>
#include <opm/input/eclipse/Schedule/Events.hpp>
@ -169,11 +170,13 @@ public:
bool isInjectionGrup(const std::string& name) const
{
OPM_TIMEFUNCTION();
return this->global_well_info.value().in_injecting_group(name);
}
bool isProductionGrup(const std::string& name) const
{
OPM_TIMEFUNCTION();
return this->global_well_info.value().in_producing_group(name);
}