OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-04 13:36:57 -06:00
Code Issues Packages Projects Releases Wiki Activity

updates

Browse Source
This commit is contained in:
Stein Krogstad 2023-10-30 20:58:09 +01:00
parent 2470b20d13
commit 7aa50f149f
11 changed files with 133 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -269,7 +269,8 @@ namespace Opm
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger,
const bool allow_switch = true) override;
const bool fixed_control = false,
const bool fixed_status = false) override;
virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
const double dt,

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

@ -1522,7 +1522,8 @@ namespace Opm
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger,
const bool allow_switch /*true*/)
const bool fixed_control /*false*/,
const bool fixed_status /*false*/)
{
//if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true;

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

@ -213,7 +213,8 @@ namespace Opm
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger,
const bool allow_switch = true) override;
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

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

@ -2321,7 +2321,8 @@ namespace Opm
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger,
const bool allow_switch /*true*/)
const bool fixed_control /*false*/,
const bool fixed_status /*false*/)
{
const int max_iter = this->param_.max_inner_iter_wells_;
@ -2339,14 +2340,21 @@ namespace Opm
const auto well_status_orig = this->wellStatus_;
auto well_status_cur = well_status_orig;
int status_switch_count = 0;
const bool allow_switching = !this->wellUnderZeroRateTarget(summary_state, well_state) && (this->well_ecl_.getStatus() == WellStatus::OPEN);
bool 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;
if (allow_switching) {
this->operability_status_.can_obtain_bhp_with_thp_limit = true;
this->operability_status_.obey_thp_limit_under_bhp_limit = true;
this->operability_status_.operable_under_only_bhp_limit = 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, allow_switch);
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++;

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

@ -535,7 +535,8 @@ intersectWithIPR(const VFPProdTable& table,
const double gfr,
const double alq,
const double ipr_a,
const double ipr_b)
const double ipr_b,
const std::function<double(const double)>& adjust_bhp)
{
// NOTE: ipr-line is q=b*bhp - a!
// ipr is given for negative flo, so
@ -547,9 +548,9 @@ intersectWithIPR(const VFPProdTable& table,
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());
auto flo_i = detail::findInterpData(ipr_a, table.getFloAxis());
detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
return std::make_pair(ipr_b, bhp_i.value);
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;
@ -558,18 +559,18 @@ intersectWithIPR(const VFPProdTable& table,
flo0 = 0.0; // start by checking flo=0
auto flo_i = detail::findInterpData(flo0, table.getFloAxis());
detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
y0 = bhp_i.value - ipr_a/ipr_b; // +0.0/ipr_b
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 = bhp_i.value + (flo1 - ipr_a)/ipr_b;
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)
w = std::clamp(w, 0.0, 1.0); // just to be safe (if y0~y1~0)
flo_x = flo0 + w*(flo1 - flo0);
}
flo0 = flo1;

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

@ -203,7 +203,8 @@ intersectWithIPR(const VFPProdTable& table,
const double gfr,
const double alq,
const double ipr_a,
const double ipr_b);
const double ipr_b,
const std::function<double(const double)>& adjust_bhp);
} // namespace detail

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

@ -90,7 +90,8 @@ double VFPProdProperties::bhp(int table_id,
const double& alq,
const double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const {
const bool use_expvfp,
const double ipr_slope) const {
const VFPProdTable& table = detail::getTable(m_tables, table_id);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg, alq, explicit_wfr,explicit_gfr, use_expvfp);
@ -167,7 +168,8 @@ EvalWell VFPProdProperties::bhp(const int table_id,
const double& alq,
const double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const
const bool use_expvfp,
const double ipr_slope /*=0*/) const
{
//Get the table
const VFPProdTable& table = detail::getTable(m_tables, table_id);
@ -192,8 +194,12 @@ EvalWell VFPProdProperties::bhp(const int table_id,
detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i);
//bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (std::max(0.0, bhp_val.dflo) * flo);
bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (bhp_val.dflo* flo);
//if (bhp_val.dflo < ipr_slope) {
// bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (std::max(0.0, bhp_val.dflo) * flo);
//} else {
bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (bhp_val.dflo* flo);
//}
bhp.setValue(bhp_val.value);
return bhp;
@ -202,7 +208,7 @@ EvalWell VFPProdProperties::bhp(const int table_id,
#define INSTANCE(...) \
template __VA_ARGS__ VFPProdProperties::bhp<__VA_ARGS__>(const int, \
const __VA_ARGS__&, const __VA_ARGS__&, const __VA_ARGS__&, \
const double&, const double&, const double&, const double&, const bool) const;
const double&, const double&, const double&, const double&, const bool, const double) const;
INSTANCE(DenseAd::Evaluation<double, -1, 4u>)
INSTANCE(DenseAd::Evaluation<double, -1, 5u>)

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

@ -68,7 +68,8 @@ public:
const double& alq,
const double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const;
const bool use_expvfp,
const double ipr_slope = 0.0) const;
/**
* Linear interpolation of bhp as a function of the input parameters
@ -90,7 +91,8 @@ public:
const double& alq,
const double& explicit_wfr,
const double& explicit_gfr,
const bool use_expvfp) const;
const bool use_expvfp,
const double ipr_slope = 0.0) const;
/**
* Linear interpolation of thp as a function of the input parameters

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

@ -365,11 +365,16 @@ 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();
double ipr_slope = 0.0;
if (use_vfpexplicit) {
const auto ipr = well_.vfpProperties()->getFloIPR(controls.vfp_table_number, well_.indexOfWell());
ipr_slope = -1/ipr.second;
}
bhp_tab = well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number,
aqua, liquid, vapour,
thp_limit,
well_.getALQ(well_state),
wfr, gfr, use_vfpexplicit);
wfr, gfr, use_vfpexplicit, ipr_slope);
}
else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + well_.name(), deferred_logger);
@ -891,12 +896,14 @@ isStableSolution(const WellState& well_state,
const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const bool use_vfpexplicit = well_.useVfpExplicit();
const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const auto bhp_adjustment = getVfpBhpAdjustment(well_state.well(well_.indexOfWell()).bhp, thp);
//const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
//const auto bhp_adjustment = getVfpBhpAdjustment(well_state.well(well_.indexOfWell()).bhp, thp);
// XXX this needs to be fixed
//assert(bhp_adjustment == 0.0);
const detail::VFPEvaluation bhp = detail::bhp(table, aqua, liquid, vapour, thp, well_.getALQ(well_state), wfr, gfr, use_vfpexplicit);
detail::VFPEvaluation bhp = detail::bhp(table, aqua, liquid, vapour, thp, well_.getALQ(well_state), wfr, gfr, use_vfpexplicit);
bhp.value = bhp.value + getVfpBhpAdjustment(bhp.value, thp);
if (bhp.dflo >= 0) {
return true;
@ -937,19 +944,20 @@ estimateStableBhp(const WellState& well_state,
auto ipr = well_.vfpProperties()->getFloIPR(controls.vfp_table_number, well_.indexOfWell());
const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const auto bhp_adjustment = getVfpBhpAdjustment(well_state.well(well_.indexOfWell()).bhp, thp);
// XXX this needs to be fixed
//assert(bhp_adjustment == 0.0);
const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth,
rho, well_.gravity());
if (ipr.first <= 0.0) {
// error message
}
const auto retval = detail::intersectWithIPR(table, thp, wfr, gfr, well_.getALQ(well_state), ipr.first+ipr.second*dp_hydro, ipr.second);
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*dp_hydro, ipr.second);
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 - dp_hydro;
return retval.value().second;
} else {
return std::nullopt;
}

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

@ -248,7 +248,8 @@ public:
const Well::ProductionControls& prod_controls,
const double WQTotal,
DeferredLogger& deferred_logger,
const bool allow_switch=true);
const bool fixed_control = false,
const bool fixed_status = false);
virtual void updatePrimaryVariables(const SummaryState& summary_state,
const WellState& well_state,
@ -416,7 +417,8 @@ protected:
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger,
const bool allow_switch = true) = 0;
const bool fixed_control = false,
const bool fixed_status = false) = 0;
virtual void updateIPRImplicit(const Simulator& ebosSimulator,
WellState& well_state,
@ -447,7 +449,12 @@ protected:
const double dt,
const double bhp,
WellState& well_state,
DeferredLogger& deferred_logger);
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);

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

@ -265,7 +265,8 @@ namespace Opm
const Well::ProductionControls& prod_controls,
const double wqTotal,
DeferredLogger& deferred_logger,
const bool allow_switch /*true*/)
const bool fixed_control,
const bool fixed_status)
{
const auto& summary_state = ebos_simulator.vanguard().summaryState();
const auto& schedule = ebos_simulator.vanguard().schedule();
@ -276,13 +277,13 @@ 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){
//std::cout << "Stopping well:" << this->name() << std::endl;
this->stopWell();
return true;
} else {
bool changed = false;
if (allow_switch) {
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);
@ -307,7 +308,7 @@ namespace Opm
}
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;
@ -334,7 +335,7 @@ namespace Opm
if (bhp_diff > 0){
//std::cout << "Re-opening well:" << this->name() << std::endl;
this->openWell();
well_state.well(this->index_of_well_).bhp = prod_limit;
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);
}
@ -342,6 +343,8 @@ namespace Opm
} else {
return false;
}
} else {
return false;
}
}
@ -499,18 +502,25 @@ namespace Opm
bool is_stable = WellBhpThpCalculator(*this).isStableSolution(well_state, this->well_ecl_, rates, summary_state, deferred_logger);
if (!is_stable) {
// solution converged to an unstable point!
//this->operability_status_.use_vfpexplicit = true;
this->operability_status_.use_vfpexplicit = true;
auto bhp_stable = WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state, deferred_logger);
if (!bhp_stable.has_value()) {
// well is not operable (no intersection)
is_operable = false;
// XXX convereged = true - iterate with closed well!
//is_operable = false;
//converged = solveWellWithZeroRate(ebos_simulator, dt, well_state, deferred_logger);
} else {
// solve equations for bhp_stable (in attempt to get closer to actual solution)
const bool converged_bhp = solveWellWithBhp(ebos_simulator, dt, bhp_stable.value(), well_state, deferred_logger);
this->updateIPRImplicit(ebos_simulator, well_state, deferred_logger);
auto bhp_stable1 = WellBhpThpCalculator(*this).estimateStableBhp(well_state, this->well_ecl_, rates, this->getRefDensity(), summary_state, deferred_logger);
const bool converged_bhp1 = solveWellWithBhp(ebos_simulator, dt, bhp_stable.value(), well_state, deferred_logger);
if (!converged_bhp) {
//
} else {
// re-solve well equations
// XXX reset thp
well_state.well(this->index_of_well_).thp = this->getTHPConstraint(summary_state);
converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
}
}
@ -520,10 +530,13 @@ namespace Opm
this->operability_status_.use_vfpexplicit = true;
//auto well_state_copy = well_state;
//
this->openWell();
auto bhp_target = estimateOperableBhp(ebos_simulator, dt, well_state, group_state, summary_state, deferred_logger);
if (!bhp_target.has_value()) {
// well can't operate using explicit fractions
is_operable = false;
// XXX convereged = true - iterate with closed well!
converged = solveWellWithZeroRate(ebos_simulator, dt, well_state, deferred_logger);
} else {
// first solve well equations for bhp = bhp_target which should give solution close to actual
converged = solveWellWithBhp(ebos_simulator, dt, bhp_target.value(), well_state, deferred_logger);
@ -531,6 +544,8 @@ namespace Opm
// debug
} else {
// re-solve well equations
// XXX reset thp
well_state.well(this->index_of_well_).thp = this->getTHPConstraint(summary_state);
converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
}
}
@ -603,12 +618,54 @@ namespace Opm
// update well-state
ws.bhp = bhp;
// solve
const bool converged = this->iterateWellEqWithSwitching(ebos_simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger, /*allow_switch*/false);
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);
/*
auto& ws = well_state.well(this->index_of_well_);
auto cmode_inj = ws.injection_cmode;
auto cmode_prod = ws.production_cmode;
if (this->isInjector()) {
assert(false);
//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::ORAT);
prod_controls.oil_rate = 0.0;
prod_controls.cmode = Well::ProducerCMode::ORAT;
ws.production_cmode = Well::ProducerCMode::ORAT;
}
*/
// 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, /*fixed_status*/ true);
this->wellStatus_ = well_status_orig;
//ws.injection_cmode = cmode_inj;
//ws.production_cmode = cmode_prod;
return converged;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::