OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-09-05 04:40:19 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into vapoilwat

Browse Source
This commit is contained in:
Paul Egberts
2022-05-11 11:55:54 +02:00
committed by GitHub
parent 6fe10a3e23 ea82bebb3a
commit a255613452
37 changed files with 753 additions and 252 deletions
Download Patch File Download Diff File Expand all files Collapse all files
opm/core/props/satfunc/RelpermDiagnostics.cpp
+13 -3
View File
@@ -112,13 +112,16 @@ namespace Opm{
const TableContainer& sof2Tables = tableManager.getSof2Tables();
const TableContainer& sgwfnTables= tableManager.getSgwfnTables();
const SwofletTable& swofletTable = tableManager.getSwofletTable();
const SgofletTable& sgofletTable = tableManager.getSgofletTable();
// Family I test.
bool family1 = pu.phase_used[BlackoilPhases::Liquid];
if (pu.phase_used[BlackoilPhases::Aqua]) {
family1 = family1 && !swofTables.empty();
family1 = family1 && (!swofTables.empty() || !swofletTable.empty());
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
family1 = family1 && (!sgofTables.empty() || !slgofTables.empty());
family1 = family1 && ((!sgofTables.empty() || !sgofletTable.empty()) || !slgofTables.empty());
}
// Family II test.
@@ -664,7 +667,9 @@ namespace Opm{
satfunc::getRawFunctionValues(tables, phases, rtep);
const TableContainer& swofTables = tables.getSwofTables();
const SwofletTable& swofletTables = tables.getSwofletTable();
const TableContainer& sgofTables = tables.getSgofTables();
const SgofletTable& sgofletTables = tables.getSgofletTable();
const TableContainer& slgofTables = tables.getSlgofTables();
const TableContainer& sof3Tables = tables.getSof3Tables();
@@ -693,14 +698,19 @@ namespace Opm{
if (!sgofTables.empty()) {
const auto& table = sgofTables.getTable<SgofTable>(satnumIdx);
krog_value = table.evaluate( "KROG" , unscaledEpsInfo_[satnumIdx].Sgl );
} else if (!sgofletTables.empty()) {
krog_value = sgofletTables[satnumIdx].krt2_relperm;
} else {
assert(!slgofTables.empty());
const auto& table = slgofTables.getTable<SlgofTable>(satnumIdx);
krog_value = table.evaluate( "KROG" , unscaledEpsInfo_[satnumIdx].Sgl );
}
{
if (!swofTables.empty()) {
const auto& table = swofTables.getTable<SwofTable>(satnumIdx);
krow_value = table.evaluate("KROW" , unscaledEpsInfo_[satnumIdx].Swl);
} else {
assert(!swofletTables.empty());
krow_value = swofletTables[satnumIdx].krt2_relperm;
}
}
if (satFamily_ == SaturationFunctionFamily::FamilyII) {
opm/simulators/flow/BlackoilModelParametersEbos.hpp
+4 -4
View File
@@ -267,7 +267,7 @@ struct MaxInnerIterWells<TypeTag, TTag::FlowModelParameters> {
};
template<class TypeTag>
struct ShutUnsolvableWells<TypeTag, TTag::FlowModelParameters> {
static constexpr bool value = false;
static constexpr bool value = true;
};
template<class TypeTag>
struct AlternativeWellRateInit<TypeTag, TTag::FlowModelParameters> {
@@ -275,7 +275,7 @@ struct AlternativeWellRateInit<TypeTag, TTag::FlowModelParameters> {
};
template<class TypeTag>
struct StrictOuterIterWells<TypeTag, TTag::FlowModelParameters> {
static constexpr int value = 99;
static constexpr int value = 6;
};
template<class TypeTag>
struct StrictInnerIterWells<TypeTag, TTag::FlowModelParameters> {
@@ -297,12 +297,12 @@ struct EnableWellOperabilityCheckIter<TypeTag, TTag::FlowModelParameters> {
template<class TypeTag>
struct RelaxedWellFlowTol<TypeTag, TTag::FlowModelParameters> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1;
static constexpr type value = 1e-3;
};
template<class TypeTag>
struct RelaxedPressureTolMsw<TypeTag, TTag::FlowModelParameters> {
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 0.5e5;
static constexpr type value = 1.0e4;
};
template<class TypeTag>
struct MaximumNumberOfWellSwitches<TypeTag, TTag::FlowModelParameters> {
opm/simulators/flow/Main.hpp
+42
View File
@@ -36,6 +36,8 @@
#include <flow/flow_ebos_brine_saltprecipitation.hpp>
#include <flow/flow_ebos_gaswater_saltprec_vapwat.hpp>
#include <flow/flow_ebos_brine_precsalt_vapwat.hpp>
#include <flow/flow_ebos_onephase.hpp>
#include <flow/flow_ebos_onephase_energy.hpp>
#include <flow/flow_ebos_oilwater_brine.hpp>
#include <flow/flow_ebos_gaswater_brine.hpp>
#include <flow/flow_ebos_energy.hpp>
@@ -291,6 +293,16 @@ private:
return this->runMICP(phases);
}
// water-only case
else if(phases.size() == 1 && phases.active(Phase::WATER) && !eclipseState_->getSimulationConfig().isThermal()) {
return this->runWaterOnly(phases);
}
// water-only case with energy
else if(phases.size() == 2 && phases.active(Phase::WATER) && eclipseState_->getSimulationConfig().isThermal()) {
return this->runWaterOnlyEnergy(phases);
}
// Twophase cases
else if (phases.size() == 2 && !eclipseState_->getSimulationConfig().isThermal()) {
return this->runTwoPhase(phases);
@@ -653,6 +665,36 @@ private:
return flowEbosFoamMain(argc_, argv_, outputCout_, outputFiles_);
}
int runWaterOnly(const Phases& phases)
{
if (!phases.active(Phase::WATER) || phases.size() != 1) {
if (outputCout_)
std::cerr << "No valid configuration is found for water-only simulation, valid options include "
<< "water, water + thermal" << std::endl;
return EXIT_FAILURE;
}
flowEbosWaterOnlySetDeck(
setupTime_, deck_, eclipseState_, schedule_, summaryConfig_);
return flowEbosWaterOnlyMain(argc_, argv_, outputCout_, outputFiles_);
}
int runWaterOnlyEnergy(const Phases& phases)
{
if (!phases.active(Phase::WATER) || phases.size() != 2) {
if (outputCout_)
std::cerr << "No valid configuration is found for water-only simulation, valid options include "
<< "water, water + thermal" << std::endl;
return EXIT_FAILURE;
}
flowEbosWaterOnlyEnergySetDeck(
setupTime_, deck_, eclipseState_, schedule_, summaryConfig_);
return flowEbosWaterOnlyEnergyMain(argc_, argv_, outputCout_, outputFiles_);
}
int runBrine(const Phases& phases)
{
if (! phases.active(Phase::WATER) || phases.size() == 2) {
opm/simulators/timestepping/ConvergenceReport.hpp
+1 -1
View File
@@ -64,7 +64,7 @@ namespace Opm
class WellFailure
{
public:
enum struct Type { Invalid, MassBalance, Pressure, ControlBHP, ControlTHP, ControlRate };
enum struct Type { Invalid, MassBalance, Pressure, ControlBHP, ControlTHP, ControlRate, Unsolvable };
WellFailure(Type t, Severity s, int phase, const std::string& well_name)
: type_(t), severity_(s), phase_(phase), well_name_(well_name)
{
opm/simulators/utils/UnsupportedFlowKeywords.cpp
+3
View File
@@ -140,8 +140,10 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
{"DPGRID", {false, std::nullopt}},
{"DPKRMOD", {false, std::nullopt}},
{"DPNUM", {false, std::nullopt}},
{"DR", {false, std::string{"Use the DRV keyword instead"}}},
{"DRILPRI", {false, std::nullopt}},
{"DSPDEINT", {false, std::nullopt}},
{"DTHETA", {false, std::string{"Use the DTHETAV keyword instead"}}},
{"DUALPERM", {false, std::nullopt}},
{"DUALPORO", {false, std::nullopt}},
{"DUMPCUPL", {false, std::nullopt}},
@@ -430,6 +432,7 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
{"OLDTRAN", {false, std::nullopt}},
{"OLDTRANR", {false, std::nullopt}},
{"OPTIONS", {false, std::nullopt}},
{"OUTRAD", {false, std::string{"Use the DRV keyword instead"}}},
{"OUTSOL", {false, std::nullopt}},
{"PARAOPTS", {false, std::nullopt}},
{"PCG32D", {false, std::nullopt}},
opm/simulators/utils/readDeck.cpp
+1 -1
View File
@@ -365,7 +365,7 @@ namespace {
if (status != MPI_SUCCESS) {
throw std::invalid_argument {
"Unable to establish cell geomtry validity across MPI ranks"
"Unable to establish cell geometry validity across MPI ranks"
};
}
#endif // HAVE_MPI
opm/simulators/wells/BlackoilWellModelGeneric.cpp
+1 -1
View File
@@ -2067,7 +2067,7 @@ forceShutWellByName(const std::string& wellname,
// process that owns it.
int well_was_shut = 0;
for (const auto& well : well_container_generic_) {
if (well->name() == wellname && !well->wellIsStopped()) {
if (well->name() == wellname) {
wellTestState().close_well(wellname, WellTestConfig::Reason::PHYSICAL, simulation_time);
well_was_shut = 1;
break;
opm/simulators/wells/BlackoilWellModel_impl.hpp
+8 -2
View File
@@ -669,7 +669,6 @@ namespace Opm {
case Well::ProducerCMode::RESV:
zero_rate_control = is_zero(prod_controls.resv_rate);
break;
default:
// Might still have zero rate controls, but is pressure controlled.
zero_rate_control = false;
@@ -1223,8 +1222,13 @@ namespace Opm {
ConvergenceReport local_report;
const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
for (const auto& well : well_container_) {
if (well->isOperableAndSolvable() ) {
if (well->isOperableAndSolvable() || well->wellIsStopped()) {
local_report += well->getWellConvergence(this->wellState(), B_avg, local_deferredLogger, iterationIdx > param_.strict_outer_iter_wells_ );
} else {
ConvergenceReport report;
using CR = ConvergenceReport;
report.setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()});
local_report += report;
}
}
@@ -1488,6 +1492,8 @@ namespace Opm {
auto& events = this->wellState().well(well->indexOfWell()).events;
if (events.hasEvent(WellState::event_mask)) {
well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), deferred_logger);
well->updatePrimaryVariables(this->wellState(), deferred_logger);
well->initPrimaryVariablesEvaluation();
// There is no new well control change input within a report step,
// so next time step, the well does not consider to have effective events anymore.
events.clearEvent(WellState::event_mask);
opm/simulators/wells/GasLiftSingleWellGeneric.cpp
+35 -12
View File
@@ -294,19 +294,38 @@ checkThpControl_() const
return thp_control;
}
std::pair<std::optional<double>,double>
GasLiftSingleWellGeneric::
computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_() const
{
auto alq = this->orig_alq_;
double new_alq = alq;
std::optional<double> bhp;
while (alq <= (this->max_alq_ + this->increment_)) {
if (bhp = computeBhpAtThpLimit_(alq); bhp) {
new_alq = alq;
break;
}
alq += this->increment_;
}
return {bhp, new_alq};
}
std::optional<GasLiftSingleWellGeneric::BasicRates>
std::pair<std::optional<GasLiftSingleWellGeneric::BasicRates>,double>
GasLiftSingleWellGeneric::
computeInitialWellRates_() const
{
std::optional<BasicRates> rates;
if (auto bhp = computeBhpAtThpLimit_(this->orig_alq_); bhp) {
double initial_alq = this->orig_alq_;
//auto alq = initial_alq;
//if (auto bhp = computeBhpAtThpLimit_(this->orig_alq_); bhp) {
if (auto [bhp, alq] = computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_(); bhp) {
{
const std::string msg = fmt::format(
"computed initial bhp {} given thp limit and given alq {}",
*bhp, this->orig_alq_);
"computed initial bhp {} given thp limit and given alq {}", *bhp, alq);
displayDebugMessage_(msg);
}
initial_alq = alq;
auto [new_bhp, bhp_is_limited] = getBhpWithLimit_(*bhp);
rates = computeWellRates_(new_bhp, bhp_is_limited);
if (rates) {
@@ -320,7 +339,7 @@ computeInitialWellRates_() const
else {
displayDebugMessage_("Aborting optimization.");
}
return rates;
return {rates, initial_alq};
}
std::optional<GasLiftSingleWellGeneric::LimitedRates>
@@ -819,12 +838,13 @@ getRateWithGroupLimit_(
}
std::optional<GasLiftSingleWellGeneric::LimitedRates>
std::pair<std::optional<GasLiftSingleWellGeneric::LimitedRates>, double>
GasLiftSingleWellGeneric::
getInitialRatesWithLimit_() const
{
std::optional<LimitedRates> limited_rates;
if (auto rates = computeInitialWellRates_(); rates) {
double initial_alq = this->orig_alq_;
if (auto [rates, alq] = computeInitialWellRates_(); rates) {
if (this->debug) {
displayDebugMessage_(
"Maybe limiting initial rates before optimize loop..");
@@ -832,8 +852,9 @@ getInitialRatesWithLimit_() const
auto temp_rates = getLimitedRatesFromRates_(*rates);
BasicRates old_rates = getWellStateRates_();
limited_rates = updateRatesToGroupLimits_(old_rates, temp_rates);
initial_alq = alq;
}
return limited_rates;
return {limited_rates, initial_alq};
}
GasLiftSingleWellGeneric::LimitedRates
@@ -1166,18 +1187,17 @@ runOptimizeLoop_(bool increase)
{
if (this->debug) debugShowProducerControlMode();
std::unique_ptr<GasLiftWellState> ret_value; // nullptr initially
auto rates = getInitialRatesWithLimit_();
auto [rates, cur_alq] = getInitialRatesWithLimit_();
if (!rates) return ret_value;
// if (this->debug) debugShowBhpAlqTable_();
if (this->debug) debugShowAlqIncreaseDecreaseCounts_();
if (this->debug) debugShowTargets_();
bool success = false; // did we succeed to increase alq?
bool alq_is_limited = false;
auto cur_alq = this->orig_alq_;
LimitedRates new_rates = *rates;
auto [temp_rates2, new_alq] = maybeAdjustALQbeforeOptimizeLoop_(
*rates, cur_alq, increase);
if (checkInitialALQmodified_(new_alq, cur_alq)) {
if (checkInitialALQmodified_(new_alq, this->orig_alq_)) {
auto delta_alq = new_alq - cur_alq;
new_rates = temp_rates2;
cur_alq = new_alq;
@@ -1591,7 +1611,10 @@ checkAlqOutsideLimits(double alq, [[maybe_unused]] double oil_rate)
// NOTE: checking for an upper limit should not be necessary
// when decreasing alq.. so this is just to catch an
// illegal state at an early point.
if (alq >= this->parent.max_alq_) {
if (this->parent.checkALQequal_(alq, this->parent.max_alq_)) {
return false;
}
else if (alq > this->parent.max_alq_) {
warn_( "unexpected: alq above upper limit when trying to "
"decrease lift gas. aborting iteration.");
result = true;
opm/simulators/wells/GasLiftSingleWellGeneric.hpp
+3 -2
View File
@@ -250,7 +250,8 @@ protected:
bool checkInitialALQmodified_(double alq, double initial_alq) const;
bool checkThpControl_() const;
virtual std::optional<double> computeBhpAtThpLimit_(double alq) const = 0;
std::optional<BasicRates> computeInitialWellRates_() const;
std::pair<std::optional<double>,double> computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_() const;
std::pair<std::optional<BasicRates>,double> computeInitialWellRates_() const;
std::optional<LimitedRates> computeLimitedWellRatesWithALQ_(double alq) const;
virtual BasicRates computeWellRates_(double bhp, bool bhp_is_limited, bool debug_output = true) const = 0;
std::optional<BasicRates> computeWellRatesWithALQ_(double alq) const;
@@ -272,7 +273,7 @@ protected:
const BasicRates& rates) const;
std::pair<double, bool> getGasRateWithGroupLimit_(
double new_gas_rate, double gas_rate) const;
std::optional<LimitedRates> getInitialRatesWithLimit_() const;
std::pair<std::optional<LimitedRates>,double> getInitialRatesWithLimit_() const;
LimitedRates getLimitedRatesFromRates_(const BasicRates& rates) const;
std::tuple<double,double,bool,bool> getLiquidRateWithGroupLimit_(
const double new_oil_rate, const double oil_rate,
opm/simulators/wells/MultisegmentWell.hpp
+3
View File
@@ -162,6 +162,9 @@ namespace Opm
protected:
int number_segments_;
// regularize msw equation
bool regularize_;
// components of the pressure drop to be included
WellSegments::CompPressureDrop compPressureDrop() const;
// multi-phase flow model
opm/simulators/wells/MultisegmentWellEval.cpp
+13 -11
View File
@@ -383,10 +383,10 @@ processFractions(const int seg) const
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
updateWellState(const BVectorWell& dwells,
const double relaxation_factor,
const double dFLimit,
const double max_pressure_change) const
updatePrimaryVariablesNewton(const BVectorWell& dwells,
const double relaxation_factor,
const double dFLimit,
const double max_pressure_change) const
{
const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_;
@@ -1622,10 +1622,10 @@ assemblePressureEq(const int seg,
}
template<typename FluidSystem, typename Indices, typename Scalar>
std::vector<Scalar>
std::pair<bool, std::vector<Scalar> >
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getWellResiduals(const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger) const
getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger) const
{
assert(int(B_avg.size() ) == baseif_.numComponents());
std::vector<Scalar> residuals(numWellEq + 1, 0.0);
@@ -1641,8 +1641,9 @@ getWellResiduals(const std::vector<Scalar>& B_avg,
}
}
if (std::isnan(residual) || std::isinf(residual)) {
OPM_DEFLOG_THROW(NumericalIssue, "nan or inf value for residal get for well " << baseif_.name()
<< " segment " << seg << " eq_idx " << eq_idx, deferred_logger);
deferred_logger.debug("nan or inf value for residal get for well " + baseif_.name()
+ " segment " + std::to_string(seg) + " eq_idx " + std::to_string(eq_idx));
return {false, residuals};
}
if (residual > residuals[eq_idx]) {
@@ -1655,12 +1656,13 @@ getWellResiduals(const std::vector<Scalar>& B_avg,
{
const double control_residual = std::abs(resWell_[0][numWellEq - 1]);
if (std::isnan(control_residual) || std::isinf(control_residual)) {
OPM_DEFLOG_THROW(NumericalIssue, "nan or inf value for control residal get for well " << baseif_.name(), deferred_logger);
deferred_logger.debug("nan or inf value for control residal get for well " + baseif_.name());
return {false, residuals};
}
residuals[numWellEq] = control_residual;
}
return residuals;
return {true, residuals};
}
template<typename FluidSystem, typename Indices, typename Scalar>
opm/simulators/wells/MultisegmentWellEval.hpp
+7 -6
View File
@@ -173,10 +173,10 @@ protected:
void updateUpwindingSegments();
// updating the well_state based on well solution dwells
void updateWellState(const BVectorWell& dwells,
const double relaxation_factor,
const double DFLimit,
const double max_pressure_change) const;
void updatePrimaryVariablesNewton(const BVectorWell& dwells,
const double relaxation_factor,
const double DFLimit,
const double max_pressure_change) const;
void computeSegmentFluidProperties(const EvalWell& temperature,
const EvalWell& saltConcentration,
@@ -199,8 +199,9 @@ protected:
EvalWell getWQTotal() const;
std::vector<Scalar> getWellResiduals(const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger) const;
std::pair<bool, std::vector<Scalar> >
getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger) const;
double getControlTolerance(const WellState& well_state,
const double tolerance_wells,
opm/simulators/wells/MultisegmentWell_impl.hpp
+32 -11
View File
@@ -49,6 +49,7 @@ namespace Opm
const std::vector<PerforationData>& perf_data)
: Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
, MSWEval(static_cast<WellInterfaceIndices<FluidSystem,Indices,Scalar>&>(*this))
, regularize_(false)
, segment_fluid_initial_(this->numberOfSegments(), std::vector<double>(this->num_components_, 0.0))
{
// not handling solvent or polymer for now with multisegment well
@@ -316,11 +317,11 @@ namespace Opm
well_potentials[phase] *= sign;
total_potential += well_potentials[phase];
}
if (total_potential < 0.0) {
if (total_potential < 0.0 && this->param_.check_well_operability_) {
// wells with negative potentials are not operable
this->operability_status_.has_negative_potentials = true;
const std::string msg = std::string("well ") + this->name() + std::string(": has non negative potentials is not operable");
deferred_logger.info(msg);
deferred_logger.warning("NEGATIVE_POTENTIALS_INOPERABLE", msg);
}
}
@@ -611,10 +612,10 @@ namespace Opm
const double dFLimit = this->param_.dwell_fraction_max_;
const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
this->MSWEval::updateWellState(dwells,
relaxation_factor,
dFLimit,
max_pressure_change);
this->MSWEval::updatePrimaryVariablesNewton(dwells,
relaxation_factor,
dFLimit,
max_pressure_change);
this->updateWellStateFromPrimaryVariables(well_state, getRefDensity(), deferred_logger);
Base::calculateReservoirRates(well_state.well(this->index_of_well_));
@@ -1373,7 +1374,14 @@ namespace Opm
const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
const WellState well_state0 = well_state;
const std::vector<Scalar> residuals0 = this->getWellResiduals(Base::B_avg_, deferred_logger);
{
// getWellFiniteResiduals returns false for nan/inf residuals
const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
if(!isFinite)
return false;
}
std::vector<std::vector<Scalar> > residual_history;
std::vector<double> measure_history;
int it = 0;
@@ -1383,14 +1391,17 @@ namespace Opm
bool converged = false;
int stagnate_count = 0;
bool relax_convergence = false;
this->regularize_ = false;
for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
if (it > this->param_.strict_inner_iter_wells_)
if (it > this->param_.strict_inner_iter_wells_) {
relax_convergence = true;
this->regularize_ = true;
}
const auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger, relax_convergence);
if (report.converged()) {
@@ -1398,12 +1409,20 @@ namespace Opm
break;
}
residual_history.push_back(this->getWellResiduals(Base::B_avg_, deferred_logger));
measure_history.push_back(this->getResidualMeasureValue(well_state,
{
// getFinteWellResiduals returns false for nan/inf residuals
const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
if (!isFinite)
return false;
residual_history.push_back(residuals);
measure_history.push_back(this->getResidualMeasureValue(well_state,
residual_history[it],
this->param_.tolerance_wells_,
this->param_.tolerance_pressure_ms_wells_,
deferred_logger) );
}
bool is_oscillate = false;
bool is_stagnate = false;
@@ -1443,6 +1462,8 @@ namespace Opm
sstr << " well " << this->name() << " observes oscillation in inner iteration " << it << "\n";
}
sstr << " relaxation_factor is " << relaxation_factor << " now\n";
this->regularize_ = true;
deferred_logger.debug(sstr.str());
}
updateWellState(dx_well, well_state, deferred_logger, relaxation_factor);
@@ -1533,7 +1554,7 @@ namespace Opm
// Add a regularization_factor to increase the accumulation term
// This will make the system less stiff and help convergence for
// difficult cases
const Scalar regularization_factor = this->param_.regularization_factor_ms_wells_;
const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_ms_wells_ : 1.0;
// for each component
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx)
opm/simulators/wells/SingleWellState.cpp
+2
View File
@@ -17,6 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/wells/SingleWellState.hpp>
#include <opm/simulators/wells/PerforationData.hpp>
@@ -39,6 +40,7 @@ SingleWellState::SingleWellState(const std::string& name_,
, surface_rates(pu_.num_phases)
, reservoir_rates(pu_.num_phases)
, perf_data(perf_input.size(), pressure_first_connection, !is_producer, pu_.num_phases)
, trivial_target(false)
{
for (std::size_t perf = 0; perf < perf_input.size(); perf++) {
this->perf_data.cell_index[perf] = perf_input[perf].cell_index;
opm/simulators/wells/SingleWellState.hpp
+1
View File
@@ -62,6 +62,7 @@ public:
std::vector<double> surface_rates;
std::vector<double> reservoir_rates;
PerfData perf_data;
bool trivial_target;
SegmentState segments;
Events events;
Well::InjectorCMode injection_cmode{Well::InjectorCMode::CMODE_UNDEFINED};
opm/simulators/wells/StandardWell_impl.hpp
+3 -2
View File
@@ -2016,11 +2016,11 @@ namespace Opm
well_potentials[phase] *= sign;
total_potential += well_potentials[phase];
}
if (total_potential < 0.0) {
if (total_potential < 0.0 && this->param_.check_well_operability_) {
// wells with negative potentials are not operable
this->operability_status_.has_negative_potentials = true;
const std::string msg = std::string("well ") + this->name() + std::string(": has negative potentials and is not operable");
deferred_logger.info(msg);
deferred_logger.warning("NEGATIVE_POTENTIALS_INOPERABLE", msg);
}
}
@@ -2476,6 +2476,7 @@ namespace Opm
// solution
std::vector<double> rates(3);
computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
this->adaptRatesForVFP(rates);
return rates;
};
opm/simulators/wells/WellInterface.hpp
+4
View File
@@ -248,6 +248,10 @@ public:
void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger);
void gliftBeginTimeStepWellTestUpdateALQ(const Simulator& ebos_simulator,
WellState& well_state,
DeferredLogger& deferred_logger);
// check whether the well is operable under the current reservoir condition
// mostly related to BHP limit and THP limit
void updateWellOperability(const Simulator& ebos_simulator,
opm/simulators/wells/WellInterfaceFluidSystem.cpp
+5 -1
View File
@@ -163,7 +163,7 @@ activeProductionConstraint(const SingleWellState& ws,
if (controls.hasControl(Well::ProducerCMode::THP) && currentControl != Well::ProducerCMode::THP) {
const auto& thp = getTHPConstraint(summaryState);
double current_thp = ws.thp;
if (thp > current_thp) {
if (thp > current_thp && !ws.trivial_target) {
// If WVFPEXP item 4 is set to YES1 or YES2
// switching to THP is prevented if the well will
// produce at a higher rate with THP control
@@ -1120,6 +1120,10 @@ getGroupProductionTargetRate(const Group& group,
const auto& rates = ws.surface_rates;
const auto current_rate = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers.
double scale = 1.0;
if (target_rate == 0.0) {
return 0.0;
}
if (current_rate > 1e-14)
scale = target_rate/current_rate;
return scale;
opm/simulators/wells/WellInterfaceGeneric.cpp
+2 -2
View File
@@ -565,9 +565,9 @@ bisectBracket(const std::function<double(const double)>& eq,
}
} else { // eq_low * eq_high > 0.0
// Still failed bracketing!
const double limit = 3.0 * unit::barsa;
const double limit = 0.1 * unit::barsa;
if (std::min(abs_low, abs_high) < limit) {
// Return the least bad solution if less off than 3 bar.
// Return the least bad solution if less off than 0.1 bar.
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE",
"Robust bhp(thp) not solved precisely for well " + this->name());
approximate_solution = abs_low < abs_high ? low : high;
opm/simulators/wells/WellInterface_impl.hpp
+56 -3
View File
@@ -332,6 +332,9 @@ namespace Opm
return;
}
if (this->isProducer()) {
gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
}
updateWellOperability(simulator, well_state_copy, deferred_logger);
if ( !this->isOperableAndSolvable() ) {
const auto msg = fmt::format("WTEST: Well {} is not operable (physical)", this->name());
@@ -446,7 +449,7 @@ namespace Opm
const GroupState& group_state,
DeferredLogger& deferred_logger)
{
if (!this->isOperableAndSolvable())
if (!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
// keep a copy of the original well state
@@ -507,6 +510,7 @@ namespace Opm
}
this->changed_to_open_this_step_ = false;
const bool well_operable = this->operability_status_.isOperableAndSolvable();
if (!well_operable && old_well_operable) {
if (this->well_ecl_.getAutomaticShutIn()) {
deferred_logger.info(" well " + this->name() + " gets SHUT during iteration ");
@@ -536,6 +540,9 @@ namespace Opm
void
WellInterface<TypeTag>::addCellRates(RateVector& rates, int cellIdx) const
{
if(!this->isOperableAndSolvable() && !this->wellIsStopped())
return;
for (int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
if (this->cells()[perfIdx] == cellIdx) {
for (int i = 0; i < RateVector::dimension; ++i) {
@@ -582,7 +589,49 @@ namespace Opm
updateWellOperability(ebos_simulator, well_state, deferred_logger);
}
template<typename TypeTag>
void
WellInterface<TypeTag>::
gliftBeginTimeStepWellTestUpdateALQ(const Simulator& ebos_simulator,
WellState& well_state,
DeferredLogger& deferred_logger)
{
const auto& summary_state = ebos_simulator.vanguard().summaryState();
const auto& well_name = this->name();
if (!this->wellHasTHPConstraints(summary_state)) {
const std::string msg = fmt::format("GLIFT WTEST: Well {} does not have THP constraints", well_name);
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& glo = schedule.glo(report_step_idx);
if (!glo.has_well(well_name)) {
const std::string msg = fmt::format(
"GLIFT WTEST: Well {} : Gas Lift not activated: "
"WLIFTOPT is probably missing. Skipping.", well_name);
deferred_logger.info(msg);
return;
}
const auto& gl_well = glo.well(well_name);
auto& max_alq_optional = gl_well.max_rate();
double max_alq;
if (max_alq_optional) {
max_alq = *max_alq_optional;
}
else {
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();
max_alq = alq_values.back();
}
well_state.setALQ(well_name, max_alq);
const std::string msg = fmt::format(
"GLIFT WTEST: Well {} : Setting ALQ to max value: {}",
well_name, max_alq);
deferred_logger.info(msg);
}
template<typename TypeTag>
void
@@ -933,6 +982,9 @@ namespace Opm
for (int p = 0; p<np; ++p) {
ws.surface_rates[p] *= scale;
}
ws.trivial_target = false;
} else {
ws.trivial_target = true;
}
break;
}
@@ -1004,8 +1056,9 @@ namespace Opm
// if more than one nonzero rate.
auto& ws = well_state.well(this->index_of_well_);
int nonzero_rate_index = -1;
const double floating_point_error_epsilon = 1e-14;
for (int p = 0; p < this->number_of_phases_; ++p) {
if (ws.surface_rates[p] != 0.0) {
if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
if (nonzero_rate_index == -1) {
nonzero_rate_index = p;
} else {