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 pull request #4188 from akva2/well_bhp_thp_calc

Browse Source 
Put well BHP/THP calculations in separate class
This commit is contained in:
Bård Skaflestad
2022-10-31 15:17:17 +01:00
committed by GitHub
parent 2840e0cc5f 796aa51e79
commit 278d7a0551
22 changed files with 1077 additions and 1171 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CMakeLists_files.cmake
+3
View File
@@ -103,6 +103,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/wells/VFPHelpers.cpp
opm/simulators/wells/VFPProdProperties.cpp
opm/simulators/wells/VFPInjProperties.cpp
opm/simulators/wells/WellBhpThpCalculator.cpp
opm/simulators/wells/WellConvergence.cpp
opm/simulators/wells/WellGroupControls.cpp
opm/simulators/wells/WellGroupHelpers.cpp
@@ -384,12 +385,14 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/wells/VFPInjProperties.hpp
opm/simulators/wells/VFPProdProperties.hpp
opm/simulators/wells/VFPProperties.hpp
opm/simulators/wells/WellBhpThpCalculator.hpp
opm/simulators/wells/WellConnectionAuxiliaryModule.hpp
opm/simulators/wells/WellConvergence.hpp
opm/simulators/wells/WellGroupControls.hpp
opm/simulators/wells/WellGroupHelpers.hpp
opm/simulators/wells/WellHelpers.hpp
opm/simulators/wells/WellInterface.hpp
opm/simulators/wells/WellInterfaceGeneric.hpp
opm/simulators/wells/WellInterface_impl.hpp
opm/simulators/wells/WellProdIndexCalculator.hpp
opm/simulators/wells/WellState.hpp
opm/simulators/wells/GasLiftSingleWell_impl.hpp
+2 -2
View File
@@ -130,8 +130,8 @@ computeBhpAtThpLimit_(double alq) const
auto bhp_at_thp_limit = this->well_.computeBhpAtThpLimitProdWithAlq(
this->ebos_simulator_,
this->summary_state_,
this->deferred_logger_,
alq);
alq,
this->deferred_logger_);
if (bhp_at_thp_limit) {
if (*bhp_at_thp_limit < this->controls_.bhp_limit) {
const std::string msg = fmt::format(
opm/simulators/wells/MultisegmentWell.hpp
+5 -5
View File
@@ -161,11 +161,11 @@ namespace Opm
double* connII,
DeferredLogger& deferred_logger) const;
virtual std::optional<double> computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
const SummaryState& summary_state,
DeferredLogger& deferred_logger,
double alq_value) const override;
std::optional<double>
computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator,
const SummaryState& summary_state,
const double alq_value,
DeferredLogger& deferred_logger) const override;
protected:
int number_segments_;
opm/simulators/wells/MultisegmentWellEval.cpp
+19 -44
View File
@@ -33,6 +33,7 @@
#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/wells/WellConvergence.hpp>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <opm/simulators/wells/WellState.hpp>
@@ -1228,7 +1229,12 @@ assembleControlEq(const WellState& well_state,
// Setup function for evaluation of BHP from THP (used only if needed).
auto bhp_from_thp = [&]() {
const auto rates = getRates();
return baseif_.calculateBhpFromThp(well_state, rates, well, summaryState, rho, deferred_logger);
return WellBhpThpCalculator(baseif_).calculateBhpFromThp(well_state,
rates,
well,
summaryState,
rho,
deferred_logger);
};
// Call generic implementation.
baseif_.assembleControlEqInj(well_state,
@@ -1246,7 +1252,12 @@ assembleControlEq(const WellState& well_state,
const auto rates = getRates();
// Setup function for evaluation of BHP from THP (used only if needed).
auto bhp_from_thp = [&]() {
return baseif_.calculateBhpFromThp(well_state, rates, well, summaryState, rho, deferred_logger);
return WellBhpThpCalculator(baseif_).calculateBhpFromThp(well_state,
rates,
well,
summaryState,
rho,
deferred_logger);
};
// Call generic implementation.
baseif_.assembleControlEqProd(well_state,
@@ -1268,47 +1279,6 @@ assembleControlEq(const WellState& well_state,
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
updateThp(WellState& well_state,
const double rho,
DeferredLogger& deferred_logger) const
{
static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua;
auto& ws = well_state.well(baseif_.indexOfWell());
// When there is no vaild VFP table provided, we set the thp to be zero.
if (!baseif_.isVFPActive(deferred_logger) || baseif_.wellIsStopped()) {
ws.thp = 0;
return;
}
// For THP controlled wells, we know the thp value
bool thp_controlled = baseif_.isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP:
ws.production_cmode == Well::ProducerCMode::THP;
if (thp_controlled) {
return;
}
// the well is under other control types, we calculate the thp based on bhp and rates
std::vector<double> rates(3, 0.0);
const PhaseUsage& pu = baseif_.phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = ws.surface_rates[pu.phase_pos[ Water ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = ws.surface_rates[pu.phase_pos[ Oil ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ];
}
ws.thp = this->calculateThpFromBhp(rates, ws.bhp, rho, deferred_logger);
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
@@ -1476,7 +1446,12 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
ws.bhp = segment_pressure[seg];
}
}
updateThp(well_state, rho, deferred_logger);
WellBhpThpCalculator(baseif_).
updateThp(rho, [this]() { return baseif_.wellEcl().alq_value(); },
{FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)},
well_state, deferred_logger);
}
template<typename FluidSystem, typename Indices, typename Scalar>
opm/simulators/wells/MultisegmentWellEval.hpp
-4
View File
@@ -221,10 +221,6 @@ protected:
// pressure drop for sub-critical valve (WSEGVALV)
EvalWell pressureDropValve(const int seg) const;
void updateThp(WellState& well_state,
const double rho,
DeferredLogger& deferred_logger) const;
void updateWellStateFromPrimaryVariables(WellState& well_state,
const double rho,
DeferredLogger& deferred_logger) const;
opm/simulators/wells/MultisegmentWellGeneric.cpp
+1 -244
View File
@@ -28,6 +28,7 @@
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/VFPHelpers.hpp>
#include <opm/simulators/wells/VFPProperties.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/wells/WellHelpers.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <opm/simulators/wells/WellState.hpp>
@@ -178,47 +179,6 @@ segmentNumberToIndex(const int segment_number) const
return segmentSet().segmentNumberToIndex(segment_number);
}
template<typename Scalar>
double
MultisegmentWellGeneric<Scalar>::
calculateThpFromBhp(const std::vector<double>& rates,
const double bhp,
const double rho,
DeferredLogger& deferred_logger) const
{
assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
const double aqua = rates[Water];
const double liquid = rates[Oil];
const double vapour = rates[Gas];
double thp = 0.0;
if (baseif_.isInjector()) {
const int table_id = baseif_.wellEcl().vfp_table_number();
const double vfp_ref_depth = baseif_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
thp = baseif_.vfpProperties()->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp);
}
else if (baseif_.isProducer()) {
const int table_id = baseif_.wellEcl().vfp_table_number();
const double alq = baseif_.wellEcl().alq_value();
const double vfp_ref_depth = baseif_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
thp = baseif_.vfpProperties()->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq);
}
else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger);
}
return thp;
}
template<typename Scalar>
void
MultisegmentWellGeneric<Scalar>::
@@ -247,209 +207,6 @@ detectOscillations(const std::vector<double>& measure_history,
stagnate = std::abs((F1 - F2) / F2) <= stagnation_rel_tol;
}
template<typename Scalar>
std::optional<double>
MultisegmentWellGeneric<Scalar>::
computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double rho,
DeferredLogger& deferred_logger) const
{
// Given a VFP function returning bhp as a function of phase
// rates and thp:
// fbhp(rates, thp),
// a function extracting the particular flow rate used for VFP
// lookups:
// flo(rates)
// and the inflow function (assuming the reservoir is fixed):
// frates(bhp)
// we want to solve the equation:
// fbhp(frates(bhp, thplimit)) - bhp = 0
// for bhp.
//
// This may result in 0, 1 or 2 solutions. If two solutions,
// the one corresponding to the lowest bhp (and therefore
// highest rate) is returned.
//
// In order to detect these situations, we will find piecewise
// linear approximations both to the inverse of the frates
// function and to the fbhp function.
//
// We first take the FLO sample points of the VFP curve, and
// find the corresponding bhp values by solving the equation:
// flo(frates(bhp)) - flo_sample = 0
// for bhp, for each flo_sample. The resulting (flo_sample,
// bhp_sample) values give a piecewise linear approximation to
// the true inverse inflow function, at the same flo values as
// the VFP data.
//
// Then we extract a piecewise linear approximation from the
// multilinear fbhp() by evaluating it at the flo_sample
// points, with fractions given by the frates(bhp_sample)
// values.
//
// When we have both piecewise linear curves defined on the
// same flo_sample points, it is easy to distinguish between
// the 0, 1 or 2 solution cases, and obtain the right interval
// in which to solve for the solution we want (with highest
// flow in case of 2 solutions).
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
// Make the fbhp() function.
const auto& controls = baseif_.wellEcl().injectionControls(summary_state);
const auto& table = baseif_.vfpProperties()->getInj()->getTable(controls.vfp_table_number);
const double vfp_ref_depth = table.getDatumDepth();
const double thp_limit = baseif_.getTHPConstraint(summary_state);
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& rates) {
assert(rates.size() == 3);
return baseif_.vfpProperties()->getInj()
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit) - dp;
};
// Make the flo() function.
auto flo = [&table](const std::vector<double>& rates) {
return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]);
};
// Get the flo samples, add extra samples at low rates and bhp
// limit point if necessary.
std::vector<double> flo_samples = table.getFloAxis();
if (flo_samples[0] > 0.0) {
const double f0 = flo_samples[0];
flo_samples.insert(flo_samples.begin(), { f0/20.0, f0/10.0, f0/5.0, f0/2.0 });
}
const double flo_bhp_limit = flo(frates(controls.bhp_limit));
if (flo_samples.back() < flo_bhp_limit) {
flo_samples.push_back(flo_bhp_limit);
}
// Find bhp values for inflow relation corresponding to flo samples.
std::vector<double> bhp_samples;
for (double flo_sample : flo_samples) {
if (flo_sample > flo_bhp_limit) {
// We would have to go over the bhp limit to obtain a
// flow of this magnitude. We associate all such flows
// with simply the bhp limit. The first one
// encountered is considered valid, the rest not. They
// are therefore skipped.
bhp_samples.push_back(controls.bhp_limit);
break;
}
auto eq = [&flo, &frates, flo_sample](double bhp) {
return flo(frates(bhp)) - flo_sample;
};
// TODO: replace hardcoded low/high limits.
const double low = 10.0 * unit::barsa;
const double high = 800.0 * unit::barsa;
const int max_iteration = 100;
const double flo_tolerance = 0.05 * std::fabs(flo_samples.back());
int iteration = 0;
try {
const double solved_bhp = RegulaFalsiBisection<WarnAndContinueOnError>::
solve(eq, low, high, max_iteration, flo_tolerance, iteration);
bhp_samples.push_back(solved_bhp);
}
catch (...) {
// Use previous value (or max value if at start) if we failed.
bhp_samples.push_back(bhp_samples.empty() ? low : bhp_samples.back());
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_EXTRACT_SAMPLES",
"Robust bhp(thp) solve failed extracting bhp values at flo samples for well " + baseif_.name());
}
}
// Find bhp values for VFP relation corresponding to flo samples.
const int num_samples = bhp_samples.size(); // Note that this can be smaller than flo_samples.size()
std::vector<double> fbhp_samples(num_samples);
for (int ii = 0; ii < num_samples; ++ii) {
fbhp_samples[ii] = fbhp(frates(bhp_samples[ii]));
}
// #define EXTRA_THP_DEBUGGING
#ifdef EXTRA_THP_DEBUGGING
std::string dbgmsg;
dbgmsg += "flo: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(flo_samples[ii]);
}
dbgmsg += "\nbhp: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(bhp_samples[ii]);
}
dbgmsg += "\nfbhp: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(fbhp_samples[ii]);
}
OpmLog::debug(dbgmsg);
#endif // EXTRA_THP_DEBUGGING
// Look for sign changes for the (fbhp_samples - bhp_samples) piecewise linear curve.
// We only look at the valid
int sign_change_index = -1;
for (int ii = 0; ii < num_samples - 1; ++ii) {
const double curr = fbhp_samples[ii] - bhp_samples[ii];
const double next = fbhp_samples[ii + 1] - bhp_samples[ii + 1];
if (curr * next < 0.0) {
// Sign change in the [ii, ii + 1] interval.
sign_change_index = ii; // May overwrite, thereby choosing the highest-flo solution.
}
}
// Handle the no solution case.
if (sign_change_index == -1) {
return std::nullopt;
}
// Solve for the proper solution in the given interval.
auto eq = [&fbhp, &frates](double bhp) {
return fbhp(frates(bhp)) - bhp;
};
// TODO: replace hardcoded low/high limits.
const double low = bhp_samples[sign_change_index + 1];
const double high = bhp_samples[sign_change_index];
const int max_iteration = 100;
const double bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0;
if (low == high) {
// We are in the high flow regime where the bhp_samples
// are all equal to the bhp_limit.
assert(low == controls.bhp_limit);
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
"Robust bhp(thp) solve failed for well " + baseif_.name());
return std::nullopt;
}
try {
const double solved_bhp = RegulaFalsiBisection<WarnAndContinueOnError>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
#ifdef EXTRA_THP_DEBUGGING
OpmLog::debug("***** " + name() + " solved_bhp = " + std::to_string(solved_bhp)
+ " flo_bhp_limit = " + std::to_string(flo_bhp_limit));
#endif // EXTRA_THP_DEBUGGING
return solved_bhp;
}
catch (...) {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
"Robust bhp(thp) solve failed for well " + baseif_.name());
return std::nullopt;
}
}
template<typename Scalar>
std::optional<double>
MultisegmentWellGeneric<Scalar>::
computeBhpAtThpLimitProdWithAlq(
const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double maxPerfPress,
const double rho,
DeferredLogger& deferred_logger,
double alq_value) const
{
return baseif_.computeBhpAtThpLimitProdCommon(frates, summary_state, maxPerfPress, rho, alq_value, deferred_logger);
}
template<typename Scalar>
bool
MultisegmentWellGeneric<Scalar>::
opm/simulators/wells/MultisegmentWellGeneric.hpp
-35
View File
@@ -60,41 +60,6 @@ protected:
/// number of segments for this well
int numberOfSegments() const;
double calculateThpFromBhp(const std::vector<double>& rates,
const double bhp,
const double rho,
DeferredLogger& deferred_logger) const;
std::optional<double> computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double rho,
DeferredLogger& deferred_logger) const;
std::optional<double> computeBhpAtThpLimitProdWithAlq(
const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double maxPerfPress,
const double rho,
DeferredLogger& deferred_logger,
double alq_value) const;
std::optional<double> bhpMax(const std::function<double(const double)>& fflo,
const double bhp_limit,
const double maxPerfPress,
const double vfp_flo_front,
DeferredLogger& deferred_logger) const;
bool bruteForceBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
DeferredLogger& deferred_logger) const;
bool bisectBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
std::optional<double>& approximate_solution,
DeferredLogger& deferred_logger) const;
/// Detect oscillation or stagnation based on the residual measure history
void detectOscillations(const std::vector<double>& measure_history,
const int it,
opm/simulators/wells/MultisegmentWell_impl.hpp
+39 -29
View File
@@ -20,6 +20,7 @@
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/input/eclipse/Schedule/MSW/Valve.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
@@ -1236,7 +1237,7 @@ namespace Opm
checkOperabilityUnderBHPLimit(const WellState& /*well_state*/, const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
{
const auto& summaryState = ebos_simulator.vanguard().summaryState();
const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
// Crude but works: default is one atmosphere.
// TODO: a better way to detect whether the BHP is defaulted or not
const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
@@ -1268,7 +1269,11 @@ namespace Opm
std::vector<double> well_rates_bhp_limit;
computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
const double thp = this->calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, getRefDensity(), deferred_logger);
const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
bhp_limit,
this->getRefDensity(),
this->wellEcl().alq_value(),
deferred_logger);
const double thp_limit = this->getTHPConstraint(summaryState);
if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
this->operability_status_.obey_thp_limit_under_bhp_limit = false;
@@ -1407,7 +1412,7 @@ namespace Opm
if (obtain_bhp) {
this->operability_status_.can_obtain_bhp_with_thp_limit = true;
const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
const double thp_limit = this->getTHPConstraint(summaryState);
@@ -1870,8 +1875,8 @@ namespace Opm
return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
ebos_simulator,
summary_state,
deferred_logger,
this->getALQ(well_state));
this->getALQ(well_state),
deferred_logger);
}
@@ -1881,8 +1886,8 @@ namespace Opm
MultisegmentWell<TypeTag>::
computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator,
const SummaryState& summary_state,
DeferredLogger& deferred_logger,
double alq_value) const
const double alq_value,
DeferredLogger& deferred_logger) const
{
// Make the frates() function.
auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
@@ -1896,15 +1901,15 @@ namespace Opm
return rates;
};
auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
computeBhpAtThpLimitProdWithAlq(frates,
auto bhpAtLimit = WellBhpThpCalculator(*this).
computeBhpAtThpLimitProd(frates,
summary_state,
maxPerfPress(ebos_simulator),
getRefDensity(),
deferred_logger,
alq_value);
this->maxPerfPress(ebos_simulator),
this->getRefDensity(),
alq_value,
deferred_logger);
if(bhpAtLimit)
if (bhpAtLimit)
return bhpAtLimit;
auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
@@ -1916,19 +1921,15 @@ namespace Opm
return rates;
};
return this->MultisegmentWellGeneric<Scalar>::
computeBhpAtThpLimitProdWithAlq(fratesIter,
return WellBhpThpCalculator(*this).
computeBhpAtThpLimitProd(fratesIter,
summary_state,
maxPerfPress(ebos_simulator),
getRefDensity(),
deferred_logger,
alq_value);
this->maxPerfPress(ebos_simulator),
this->getRefDensity(),
alq_value,
deferred_logger);
}
template<typename TypeTag>
std::optional<double>
MultisegmentWell<TypeTag>::
@@ -1948,13 +1949,16 @@ namespace Opm
return rates;
};
auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
auto bhpAtLimit = WellBhpThpCalculator(*this).
computeBhpAtThpLimitInj(frates,
summary_state,
getRefDensity(),
this->getRefDensity(),
0.05,
100,
false,
deferred_logger);
if(bhpAtLimit)
if (bhpAtLimit)
return bhpAtLimit;
auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
@@ -1966,8 +1970,14 @@ namespace Opm
return rates;
};
return this->MultisegmentWellGeneric<Scalar>::
computeBhpAtThpLimitInj(fratesIter, summary_state, getRefDensity(), deferred_logger);
return WellBhpThpCalculator(*this).
computeBhpAtThpLimitInj(fratesIter,
summary_state,
this->getRefDensity(),
0.05,
100,
false,
deferred_logger);
}
opm/simulators/wells/StandardWell.hpp
+3 -3
View File
@@ -218,11 +218,11 @@ namespace Opm
std::vector<double> &potentials,
double alq) const;
virtual std::optional<double> computeBhpAtThpLimitProdWithAlq(
std::optional<double> computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
const SummaryState& summary_state,
DeferredLogger& deferred_logger,
double alq_value) const override;
const double alq_value,
DeferredLogger& deferred_logger) const override;
virtual void computeWellRatesWithBhp(
const Simulator& ebosSimulator,
opm/simulators/wells/StandardWellEval.cpp
+20 -48
View File
@@ -32,6 +32,7 @@
#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/ParallelWellInfo.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/wells/WellConvergence.hpp>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <opm/simulators/wells/WellState.hpp>
@@ -393,7 +394,12 @@ assembleControlEq(const WellState& well_state,
// Setup function for evaluation of BHP from THP (used only if needed).
auto bhp_from_thp = [&]() {
const auto rates = getRates();
return baseif_.calculateBhpFromThp(well_state, rates, well, summaryState, this->getRho(), deferred_logger);
return WellBhpThpCalculator(baseif_).calculateBhpFromThp(well_state,
rates,
well,
summaryState,
this->getRho(),
deferred_logger);
};
// Call generic implementation.
const auto& inj_controls = well.injectionControls(summaryState);
@@ -412,7 +418,12 @@ assembleControlEq(const WellState& well_state,
const auto rates = getRates();
// Setup function for evaluation of BHP from THP (used only if needed).
auto bhp_from_thp = [&]() {
return baseif_.calculateBhpFromThp(well_state, rates, well, summaryState, this->getRho(), deferred_logger);
return WellBhpThpCalculator(baseif_).calculateBhpFromThp(well_state,
rates,
well,
summaryState,
this->getRho(),
deferred_logger);
};
// Call generic implementation.
const auto& prod_controls = well.productionControls(summaryState);
@@ -556,51 +567,6 @@ processFractions() const
}
}
template<class FluidSystem, class Indices, class Scalar>
void
StandardWellEval<FluidSystem,Indices,Scalar>::
updateThp(WellState& well_state,
DeferredLogger& deferred_logger) const
{
static constexpr int Gas = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Gas;
static constexpr int Oil = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Oil;
static constexpr int Water = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Water;
auto& ws = well_state.well(baseif_.indexOfWell());
// When there is no vaild VFP table provided, we set the thp to be zero.
if (!baseif_.isVFPActive(deferred_logger) || baseif_.wellIsStopped()) {
ws.thp = 0;
return;
}
// For THP controlled wells, we know the thp value
bool thp_controlled = baseif_.isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP:
ws.production_cmode == Well::ProducerCMode::THP;
if (thp_controlled) {
return;
}
// the well is under other control types, we calculate the thp based on bhp and rates
std::vector<double> rates(3, 0.0);
const PhaseUsage& pu = baseif_.phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = ws.surface_rates[pu.phase_pos[ Water ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = ws.surface_rates[pu.phase_pos[ Oil ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ];
}
ws.thp = this->calculateThpFromBhp(well_state,
rates,
ws.bhp,
deferred_logger);
}
template<class FluidSystem, class Indices, class Scalar>
void
StandardWellEval<FluidSystem,Indices,Scalar>::
@@ -702,7 +668,13 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
}
}
updateThp(well_state, deferred_logger);
WellBhpThpCalculator(baseif_).
updateThp(this->getRho(),
[this,&well_state]() { return this->baseif_.getALQ(well_state); },
{FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx),
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)},
well_state, deferred_logger);
}
template<class FluidSystem, class Indices, class Scalar>
opm/simulators/wells/StandardWellEval.hpp
-3
View File
@@ -177,9 +177,6 @@ protected:
void updateWellStateFromPrimaryVariablesPolyMW(WellState& well_state) const;
void updateThp(WellState& well_state,
DeferredLogger& deferred_logger) const;
// total number of the well equations and primary variables
// there might be extra equations be used, numWellEq will be updated during the initialization
int numWellEq_ = numStaticWellEq;
opm/simulators/wells/StandardWellGeneric.cpp
+1 -242
View File
@@ -35,6 +35,7 @@
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/VFPHelpers.hpp>
#include <opm/simulators/wells/VFPProperties.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/wells/WellHelpers.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <opm/simulators/wells/WellState.hpp>
@@ -112,48 +113,6 @@ relaxationFactorFraction(const double old_value,
return relaxation_factor;
}
template<class Scalar>
double
StandardWellGeneric<Scalar>::
calculateThpFromBhp(const WellState &well_state,
const std::vector<double>& rates,
const double bhp,
DeferredLogger& deferred_logger) const
{
assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
const double aqua = rates[Water];
const double liquid = rates[Oil];
const double vapour = rates[Gas];
// pick the density in the top layer
double thp = 0.0;
if (baseif_.isInjector()) {
const int table_id = baseif_.wellEcl().vfp_table_number();
const double vfp_ref_depth = baseif_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, getRho(), baseif_.gravity());
thp = baseif_.vfpProperties()->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp);
}
else if (baseif_.isProducer()) {
const int table_id = baseif_.wellEcl().vfp_table_number();
const double alq = baseif_.getALQ(well_state);
const double vfp_ref_depth = baseif_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, getRho(), baseif_.gravity());
thp = baseif_.vfpProperties()->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq);
}
else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger);
}
return thp;
}
template<class Scalar>
void
StandardWellGeneric<Scalar>::
@@ -191,206 +150,6 @@ computeConnectionPressureDelta()
baseif_.parallelWellInfo().partialSumPerfValues(beg, end);
}
template<class Scalar>
std::optional<double>
StandardWellGeneric<Scalar>::
computeBhpAtThpLimitProdWithAlq(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
DeferredLogger& deferred_logger,
double maxPerfPress,
double alq_value) const
{
return baseif_.computeBhpAtThpLimitProdCommon(frates, summary_state, maxPerfPress, getRho(), alq_value, deferred_logger);
}
template<class Scalar>
std::optional<double>
StandardWellGeneric<Scalar>::
computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const
{
// Given a VFP function returning bhp as a function of phase
// rates and thp:
// fbhp(rates, thp),
// a function extracting the particular flow rate used for VFP
// lookups:
// flo(rates)
// and the inflow function (assuming the reservoir is fixed):
// frates(bhp)
// we want to solve the equation:
// fbhp(frates(bhp, thplimit)) - bhp = 0
// for bhp.
//
// This may result in 0, 1 or 2 solutions. If two solutions,
// the one corresponding to the lowest bhp (and therefore
// highest rate) is returned.
//
// In order to detect these situations, we will find piecewise
// linear approximations both to the inverse of the frates
// function and to the fbhp function.
//
// We first take the FLO sample points of the VFP curve, and
// find the corresponding bhp values by solving the equation:
// flo(frates(bhp)) - flo_sample = 0
// for bhp, for each flo_sample. The resulting (flo_sample,
// bhp_sample) values give a piecewise linear approximation to
// the true inverse inflow function, at the same flo values as
// the VFP data.
//
// Then we extract a piecewise linear approximation from the
// multilinear fbhp() by evaluating it at the flo_sample
// points, with fractions given by the frates(bhp_sample)
// values.
//
// When we have both piecewise linear curves defined on the
// same flo_sample points, it is easy to distinguish between
// the 0, 1 or 2 solution cases, and obtain the right interval
// in which to solve for the solution we want (with highest
// flow in case of 2 solutions).
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
// Make the fbhp() function.
const auto& controls = baseif_.wellEcl().injectionControls(summary_state);
const auto& table = baseif_.vfpProperties()->getInj()->getTable(controls.vfp_table_number);
const double vfp_ref_depth = table.getDatumDepth();
const double thp_limit = baseif_.getTHPConstraint(summary_state);
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, getRho(), baseif_.gravity());
auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& rates) {
assert(rates.size() == 3);
return baseif_.vfpProperties()->getInj()
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit) - dp;
};
// Make the flo() function.
auto flo = [&table](const std::vector<double>& rates) {
return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]);
};
// Get the flo samples, add extra samples at low rates and bhp
// limit point if necessary.
std::vector<double> flo_samples = table.getFloAxis();
if (flo_samples[0] > 0.0) {
const double f0 = flo_samples[0];
flo_samples.insert(flo_samples.begin(), { f0/20.0, f0/10.0, f0/5.0, f0/2.0 });
}
const double flo_bhp_limit = flo(frates(controls.bhp_limit));
if (flo_samples.back() < flo_bhp_limit) {
flo_samples.push_back(flo_bhp_limit);
}
// Find bhp values for inflow relation corresponding to flo samples.
std::vector<double> bhp_samples;
for (double flo_sample : flo_samples) {
if (flo_sample > flo_bhp_limit) {
// We would have to go over the bhp limit to obtain a
// flow of this magnitude. We associate all such flows
// with simply the bhp limit. The first one
// encountered is considered valid, the rest not. They
// are therefore skipped.
bhp_samples.push_back(controls.bhp_limit);
break;
}
auto eq = [&flo, &frates, flo_sample](double bhp) {
return flo(frates(bhp)) - flo_sample;
};
// TODO: replace hardcoded low/high limits.
const double low = 10.0 * unit::barsa;
const double high = 800.0 * unit::barsa;
const int max_iteration = 50;
const double flo_tolerance = 1e-6 * std::fabs(flo_samples.back());
int iteration = 0;
try {
const double solved_bhp = RegulaFalsiBisection<>::
solve(eq, low, high, max_iteration, flo_tolerance, iteration);
bhp_samples.push_back(solved_bhp);
}
catch (...) {
// Use previous value (or max value if at start) if we failed.
bhp_samples.push_back(bhp_samples.empty() ? low : bhp_samples.back());
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_EXTRACT_SAMPLES",
"Robust bhp(thp) solve failed extracting bhp values at flo samples for well " + baseif_.name());
}
}
// Find bhp values for VFP relation corresponding to flo samples.
const int num_samples = bhp_samples.size(); // Note that this can be smaller than flo_samples.size()
std::vector<double> fbhp_samples(num_samples);
for (int ii = 0; ii < num_samples; ++ii) {
fbhp_samples[ii] = fbhp(frates(bhp_samples[ii]));
}
// #define EXTRA_THP_DEBUGGING
#ifdef EXTRA_THP_DEBUGGING
std::string dbgmsg;
dbgmsg += "flo: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(flo_samples[ii]);
}
dbgmsg += "\nbhp: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(bhp_samples[ii]);
}
dbgmsg += "\nfbhp: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(fbhp_samples[ii]);
}
OpmLog::debug(dbgmsg);
#endif // EXTRA_THP_DEBUGGING
// Look for sign changes for the (fbhp_samples - bhp_samples) piecewise linear curve.
// We only look at the valid
int sign_change_index = -1;
for (int ii = 0; ii < num_samples - 1; ++ii) {
const double curr = fbhp_samples[ii] - bhp_samples[ii];
const double next = fbhp_samples[ii + 1] - bhp_samples[ii + 1];
if (curr * next < 0.0) {
// Sign change in the [ii, ii + 1] interval.
sign_change_index = ii; // May overwrite, thereby choosing the highest-flo solution.
}
}
// Handle the no solution case.
if (sign_change_index == -1) {
return std::nullopt;
}
// Solve for the proper solution in the given interval.
auto eq = [&fbhp, &frates](double bhp) {
return fbhp(frates(bhp)) - bhp;
};
// TODO: replace hardcoded low/high limits.
const double low = bhp_samples[sign_change_index + 1];
const double high = bhp_samples[sign_change_index];
const int max_iteration = 50;
const double bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0;
if (low == high) {
// We are in the high flow regime where the bhp_samples
// are all equal to the bhp_limit.
assert(low == controls.bhp_limit);
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
"Robust bhp(thp) solve failed for well " + baseif_.name());
return std::nullopt;
}
try {
const double solved_bhp = RegulaFalsiBisection<>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
#ifdef EXTRA_THP_DEBUGGING
OpmLog::debug("***** " + name() + " solved_bhp = " + std::to_string(solved_bhp)
+ " flo_bhp_limit = " + std::to_string(flo_bhp_limit));
#endif // EXTRA_THP_DEBUGGING
return solved_bhp;
}
catch (...) {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
"Robust bhp(thp) solve failed for well " + baseif_.name());
return std::nullopt;
}
}
template<class Scalar>
unsigned int
StandardWellGeneric<Scalar>::
opm/simulators/wells/StandardWellGeneric.hpp
+4 -15
View File
@@ -79,22 +79,8 @@ protected:
static double relaxationFactorFraction(const double old_value,
const double dx);
double calculateThpFromBhp(const WellState& well_state,
const std::vector<double>& rates,
const double bhp,
DeferredLogger& deferred_logger) const;
void computeConnectionPressureDelta();
std::optional<double> computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const;
std::optional<double> computeBhpAtThpLimitProdWithAlq(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
DeferredLogger& deferred_logger,
double maxPerfPress,
double alq_value) const;
// Base interface reference
const WellInterfaceGeneric& baseif_;
@@ -120,7 +106,10 @@ protected:
mutable BVectorWell Bx_;
mutable BVectorWell invDrw_;
double getRho() const { return perf_densities_[0]; }
double getRho() const
{
return this->perf_densities_.empty() ? 0.0 : perf_densities_[0];
}
};
}
opm/simulators/wells/StandardWell_impl.hpp
+34 -23
View File
@@ -23,6 +23,7 @@
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/linalg/SmallDenseMatrixUtils.hpp>
#include <opm/simulators/wells/VFPHelpers.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/wells/WellConvergence.hpp>
#include <algorithm>
@@ -1102,7 +1103,7 @@ namespace Opm
checkOperabilityUnderBHPLimit(const WellState& well_state, const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
{
const auto& summaryState = ebos_simulator.vanguard().summaryState();
const double bhp_limit = this->mostStrictBhpFromBhpLimits(summaryState);
const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
// Crude but works: default is one atmosphere.
// TODO: a better way to detect whether the BHP is defaulted or not
const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
@@ -1135,7 +1136,11 @@ namespace Opm
computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
this->adaptRatesForVFP(well_rates_bhp_limit);
const double thp = this->calculateThpFromBhp(well_state, well_rates_bhp_limit, bhp_limit, deferred_logger);
const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
bhp_limit,
this->getRho(),
this->getALQ(well_state),
deferred_logger);
const double thp_limit = this->getTHPConstraint(summaryState);
if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
this->operability_status_.obey_thp_limit_under_bhp_limit = false;
@@ -1170,7 +1175,7 @@ namespace Opm
if (obtain_bhp) {
this->operability_status_.can_obtain_bhp_with_thp_limit = true;
const double bhp_limit = this->mostStrictBhpFromBhpLimits(summaryState);
const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
const double thp_limit = this->getTHPConstraint(summaryState);
@@ -1924,7 +1929,7 @@ namespace Opm
{
double bhp;
auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
ebos_simulator, summary_state, deferred_logger, alq);
ebos_simulator, summary_state, alq, deferred_logger);
if (bhp_at_thp_limit) {
const auto& controls = this->well_ecl_.productionControls(summary_state);
bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
@@ -2017,7 +2022,7 @@ namespace Opm
const auto& summaryState = ebosSimulator.vanguard().summaryState();
if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
// get the bhp value based on the bhp constraints
double bhp = this->mostStrictBhpFromBhpLimits(summaryState);
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
@@ -2564,8 +2569,8 @@ namespace Opm
{
return computeBhpAtThpLimitProdWithAlq(ebos_simulator,
summary_state,
deferred_logger,
this->getALQ(well_state));
this->getALQ(well_state),
deferred_logger);
}
template<typename TypeTag>
@@ -2573,8 +2578,8 @@ namespace Opm
StandardWell<TypeTag>::
computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator,
const SummaryState& summary_state,
DeferredLogger& deferred_logger,
double alq_value) const
const double alq_value,
DeferredLogger& deferred_logger) const
{
// Make the frates() function.
auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
@@ -2597,13 +2602,14 @@ namespace Opm
double pressure_cell = this->getPerfCellPressure(fs).value();
max_pressure = std::max(max_pressure, pressure_cell);
}
auto bhpAtLimit = this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitProdWithAlq(frates,
summary_state,
deferred_logger,
max_pressure,
alq_value);
auto bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(frates,
summary_state,
max_pressure,
this->getRho(),
alq_value,
deferred_logger);
auto v = frates(*bhpAtLimit);
if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }))
if (bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }))
return bhpAtLimit;
auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
@@ -2616,11 +2622,12 @@ namespace Opm
return rates;
};
bhpAtLimit = this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitProdWithAlq(fratesIter,
summary_state,
deferred_logger,
max_pressure,
alq_value);
bhpAtLimit = WellBhpThpCalculator(*this).computeBhpAtThpLimitProd(fratesIter,
summary_state,
max_pressure,
this->getRho(),
alq_value,
deferred_logger);
v = frates(*bhpAtLimit);
if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }))
return bhpAtLimit;
@@ -2650,9 +2657,13 @@ namespace Opm
return rates;
};
return this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitInj(frates,
summary_state,
deferred_logger);
return WellBhpThpCalculator(*this).computeBhpAtThpLimitInj(frates,
summary_state,
this->getRho(),
1e-6,
50,
true,
deferred_logger);
}
opm/simulators/wells/WellBhpThpCalculator.cpp
+786
View File
@@ -0,0 +1,786 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2018 IRIS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/common/utility/numeric/RootFinders.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/input/eclipse/Schedule/VFPInjTable.hpp>
#include <opm/input/eclipse/Schedule/Well/Well.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/VFPProperties.hpp>
#include <opm/simulators/wells/WellHelpers.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <cassert>
#include <cmath>
static constexpr bool extraBhpAtThpLimitOutput = false;
namespace Opm
{
bool
WellBhpThpCalculator::wellHasTHPConstraints(const SummaryState& summaryState) const
{
const auto& well_ecl = well_.wellEcl();
if (well_ecl.isInjector()) {
const auto controls = well_ecl.injectionControls(summaryState);
if (controls.hasControl(Well::InjectorCMode::THP))
return true;
}
if (well_ecl.isProducer()) {
const auto controls = well_ecl.productionControls(summaryState);
if (controls.hasControl(Well::ProducerCMode::THP))
return true;
}
return false;
}
double WellBhpThpCalculator::getTHPConstraint(const SummaryState& summaryState) const
{
const auto& well_ecl = well_.wellEcl();
if (well_ecl.isInjector()) {
const auto& controls = well_ecl.injectionControls(summaryState);
return controls.thp_limit;
}
if (well_ecl.isProducer( )) {
const auto& controls = well_ecl.productionControls(summaryState);
return controls.thp_limit;
}
return 0.0;
}
double WellBhpThpCalculator::mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const
{
const auto& well_ecl = well_.wellEcl();
if (well_ecl.isInjector()) {
const auto& controls = well_ecl.injectionControls(summaryState);
return controls.bhp_limit;
}
if (well_ecl.isProducer( )) {
const auto& controls = well_ecl.productionControls(summaryState);
return controls.bhp_limit;
}
return 0.0;
}
double WellBhpThpCalculator::calculateThpFromBhp(const std::vector<double>& rates,
const double bhp,
const double rho,
const double alq,
DeferredLogger& deferred_logger) const
{
assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
const double aqua = rates[Water];
const double liquid = rates[Oil];
const double vapour = rates[Gas];
// pick the density in the top layer
double thp = 0.0;
if (well_.isInjector()) {
const int table_id = well_.wellEcl().vfp_table_number();
const double vfp_ref_depth = well_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
thp = well_.vfpProperties()->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp);
}
else if (well_.isProducer()) {
const int table_id = well_.wellEcl().vfp_table_number();
const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
thp = well_.vfpProperties()->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq);
}
else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger);
}
return thp;
}
std::optional<double>
WellBhpThpCalculator::
computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double maxPerfPress,
const double rho,
const double alq_value,
DeferredLogger& deferred_logger) const
{
// Given a VFP function returning bhp as a function of phase
// rates and thp:
// fbhp(rates, thp),
// a function extracting the particular flow rate used for VFP
// lookups:
// flo(rates)
// and the inflow function (assuming the reservoir is fixed):
// frates(bhp)
// we want to solve the equation:
// fbhp(frates(bhp, thplimit)) - bhp = 0
// for bhp.
//
// This may result in 0, 1 or 2 solutions. If two solutions,
// the one corresponding to the lowest bhp (and therefore
// highest rate) should be returned.
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
// Make the fbhp() function.
const auto& controls = well_.wellEcl().productionControls(summary_state);
const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const double vfp_ref_depth = table.getDatumDepth();
const double thp_limit = this->getTHPConstraint(summary_state);
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
auto fbhp = [this, &controls, thp_limit, dp, alq_value](const std::vector<double>& rates) {
assert(rates.size() == 3);
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_vfpexp = well_.useVfpExplicit();
return well_.vfpProperties()->getProd()
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, alq_value, wfr, gfr, use_vfpexp) - dp;
};
// Make the flo() function.
auto flo = [&table](const std::vector<double>& rates) {
return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]);
};
// Find the bhp-point where production becomes nonzero.
auto fflo = [&flo, &frates](double bhp) { return flo(frates(bhp)); };
auto bhp_max = this->bhpMax(fflo, controls.bhp_limit, maxPerfPress, table.getFloAxis().front(), deferred_logger);
// could not solve for the bhp-point, we could not continue to find the bhp
if (!bhp_max.has_value()) {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE",
"Robust bhp(thp) solve failed due to not being able to "
"find bhp-point where production becomes non-zero for well " + well_.name());
return std::nullopt;
}
const std::array<double, 2> range {controls.bhp_limit, *bhp_max};
return this->computeBhpAtThpLimit(frates, fbhp, range, deferred_logger);
}
std::optional<double>
WellBhpThpCalculator::
computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double rho,
const double flo_rel_tol,
const int max_iteration,
const bool throwOnError,
DeferredLogger& deferred_logger) const
{
if (throwOnError) {
return computeBhpAtThpLimitInjImpl<ThrowOnError>(frates, summary_state,
rho, flo_rel_tol,
max_iteration, deferred_logger);
} else {
return computeBhpAtThpLimitInjImpl<WarnAndContinueOnError>(frates, summary_state,
rho, flo_rel_tol,
max_iteration, deferred_logger);
}
}
void WellBhpThpCalculator::updateThp(const double rho,
const std::function<double()>& alq_value,
const std::array<unsigned,3>& active,
WellState& well_state,
DeferredLogger& deferred_logger) const
{
static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua;
auto& ws = well_state.well(well_.indexOfWell());
// When there is no vaild VFP table provided, we set the thp to be zero.
if (!well_.isVFPActive(deferred_logger) || well_.wellIsStopped()) {
ws.thp = 0;
return;
}
// For THP controlled wells, we know the thp value
bool thp_controlled = well_.isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP:
ws.production_cmode == Well::ProducerCMode::THP;
if (thp_controlled) {
return;
}
// the well is under other control types, we calculate the thp based on bhp and rates
std::vector<double> rates(3, 0.0);
const PhaseUsage& pu = well_.phaseUsage();
if (active[Water]) {
rates[ Water ] = ws.surface_rates[pu.phase_pos[ Water ] ];
}
if (active[Oil]) {
rates[ Oil ] = ws.surface_rates[pu.phase_pos[ Oil ] ];
}
if (active[Gas]) {
rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ];
}
ws.thp = this->calculateThpFromBhp(rates, ws.bhp, rho, alq_value(), deferred_logger);
}
template<class EvalWell>
EvalWell WellBhpThpCalculator::
calculateBhpFromThp(const WellState& well_state,
const std::vector<EvalWell>& rates,
const Well& well,
const SummaryState& summaryState,
const double rho,
DeferredLogger& deferred_logger) const
{
// TODO: when well is under THP control, the BHP is dependent on the rates,
// the well rates is also dependent on the BHP, so it might need to do some iteration.
// However, when group control is involved, change of the rates might impacts other wells
// so iterations on a higher level will be required. Some investigation might be needed when
// we face problems under THP control.
assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua;
const EvalWell aqua = rates[Water];
const EvalWell liquid = rates[Oil];
const EvalWell vapour = rates[Gas];
// pick the reference density
// typically the reference in the top layer
if (well_.isInjector() )
{
const auto& controls = well.injectionControls(summaryState);
const double vfp_ref_depth = well_.vfpProperties()->getInj()->getTable(controls.vfp_table_number).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
return well_.vfpProperties()->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, well_.getTHPConstraint(summaryState)) - dp;
}
else if (well_.isProducer()) {
const auto& controls = well.productionControls(summaryState);
const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
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();
return well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number,
aqua, liquid, vapour,
well_.getTHPConstraint(summaryState),
well_.getALQ(well_state),
wfr, gfr, use_vfpexplicit) - dp;
}
else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + well_.name(), deferred_logger);
}
}
template<class ErrorPolicy>
std::optional<double>
WellBhpThpCalculator::
computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double rho,
const double flo_rel_tol,
const int max_iteration,
DeferredLogger& deferred_logger) const
{
// Given a VFP function returning bhp as a function of phase
// rates and thp:
// fbhp(rates, thp),
// a function extracting the particular flow rate used for VFP
// lookups:
// flo(rates)
// and the inflow function (assuming the reservoir is fixed):
// frates(bhp)
// we want to solve the equation:
// fbhp(frates(bhp, thplimit)) - bhp = 0
// for bhp.
//
// This may result in 0, 1 or 2 solutions. If two solutions,
// the one corresponding to the lowest bhp (and therefore
// highest rate) is returned.
//
// In order to detect these situations, we will find piecewise
// linear approximations both to the inverse of the frates
// function and to the fbhp function.
//
// We first take the FLO sample points of the VFP curve, and
// find the corresponding bhp values by solving the equation:
// flo(frates(bhp)) - flo_sample = 0
// for bhp, for each flo_sample. The resulting (flo_sample,
// bhp_sample) values give a piecewise linear approximation to
// the true inverse inflow function, at the same flo values as
// the VFP data.
//
// Then we extract a piecewise linear approximation from the
// multilinear fbhp() by evaluating it at the flo_sample
// points, with fractions given by the frates(bhp_sample)
// values.
//
// When we have both piecewise linear curves defined on the
// same flo_sample points, it is easy to distinguish between
// the 0, 1 or 2 solution cases, and obtain the right interval
// in which to solve for the solution we want (with highest
// flow in case of 2 solutions).
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
// Make the fbhp() function.
const auto& controls = well_.wellEcl().injectionControls(summary_state);
const auto& table = well_.vfpProperties()->getInj()->getTable(controls.vfp_table_number);
const double vfp_ref_depth = table.getDatumDepth();
const double thp_limit = this->getTHPConstraint(summary_state);
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(),
vfp_ref_depth,
rho, well_.gravity());
auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& rates) {
assert(rates.size() == 3);
return well_.vfpProperties()->getInj()
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit) - dp;
};
// Make the flo() function.
auto flo = [&table](const std::vector<double>& rates) {
return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]);
};
// Get the flo samples, add extra samples at low rates and bhp
// limit point if necessary.
std::vector<double> flo_samples = table.getFloAxis();
if (flo_samples[0] > 0.0) {
const double f0 = flo_samples[0];
flo_samples.insert(flo_samples.begin(), { f0/20.0, f0/10.0, f0/5.0, f0/2.0 });
}
const double flo_bhp_limit = flo(frates(controls.bhp_limit));
if (flo_samples.back() < flo_bhp_limit) {
flo_samples.push_back(flo_bhp_limit);
}
// Find bhp values for inflow relation corresponding to flo samples.
std::vector<double> bhp_samples;
for (double flo_sample : flo_samples) {
if (flo_sample > flo_bhp_limit) {
// We would have to go over the bhp limit to obtain a
// flow of this magnitude. We associate all such flows
// with simply the bhp limit. The first one
// encountered is considered valid, the rest not. They
// are therefore skipped.
bhp_samples.push_back(controls.bhp_limit);
break;
}
auto eq = [&flo, &frates, flo_sample](double bhp) {
return flo(frates(bhp)) - flo_sample;
};
// TODO: replace hardcoded low/high limits.
const double low = 10.0 * unit::barsa;
const double high = 800.0 * unit::barsa;
const double flo_tolerance = flo_rel_tol * std::fabs(flo_samples.back());
int iteration = 0;
try {
const double solved_bhp = RegulaFalsiBisection<ErrorPolicy>::
solve(eq, low, high, max_iteration, flo_tolerance, iteration);
bhp_samples.push_back(solved_bhp);
}
catch (...) {
// Use previous value (or max value if at start) if we failed.
bhp_samples.push_back(bhp_samples.empty() ? low : bhp_samples.back());
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_EXTRACT_SAMPLES",
"Robust bhp(thp) solve failed extracting bhp values at flo samples for well " + well_.name());
}
}
// Find bhp values for VFP relation corresponding to flo samples.
const int num_samples = bhp_samples.size(); // Note that this can be smaller than flo_samples.size()
std::vector<double> fbhp_samples(num_samples);
for (int ii = 0; ii < num_samples; ++ii) {
fbhp_samples[ii] = fbhp(frates(bhp_samples[ii]));
}
if constexpr (extraBhpAtThpLimitOutput) {
std::string dbgmsg;
dbgmsg += "flo: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(flo_samples[ii]);
}
dbgmsg += "\nbhp: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(bhp_samples[ii]);
}
dbgmsg += "\nfbhp: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(fbhp_samples[ii]);
}
OpmLog::debug(dbgmsg);
}
// Look for sign changes for the (fbhp_samples - bhp_samples) piecewise linear curve.
// We only look at the valid
int sign_change_index = -1;
for (int ii = 0; ii < num_samples - 1; ++ii) {
const double curr = fbhp_samples[ii] - bhp_samples[ii];
const double next = fbhp_samples[ii + 1] - bhp_samples[ii + 1];
if (curr * next < 0.0) {
// Sign change in the [ii, ii + 1] interval.
sign_change_index = ii; // May overwrite, thereby choosing the highest-flo solution.
}
}
// Handle the no solution case.
if (sign_change_index == -1) {
return std::nullopt;
}
// Solve for the proper solution in the given interval.
auto eq = [&fbhp, &frates](double bhp) {
return fbhp(frates(bhp)) - bhp;
};
// TODO: replace hardcoded low/high limits.
const double low = bhp_samples[sign_change_index + 1];
const double high = bhp_samples[sign_change_index];
const double bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0;
if (low == high) {
// We are in the high flow regime where the bhp_samples
// are all equal to the bhp_limit.
assert(low == controls.bhp_limit);
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
"Robust bhp(thp) solve failed for well " + well_.name());
return std::nullopt;
}
try {
const double solved_bhp = RegulaFalsiBisection<ErrorPolicy>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
if constexpr (extraBhpAtThpLimitOutput) {
OpmLog::debug("***** " + well_.name() + " solved_bhp = " + std::to_string(solved_bhp)
+ " flo_bhp_limit = " + std::to_string(flo_bhp_limit));
}
return solved_bhp;
}
catch (...) {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
"Robust bhp(thp) solve failed for well " + well_.name());
return std::nullopt;
}
}
std::optional<double>
WellBhpThpCalculator::
bhpMax(const std::function<double(const double)>& fflo,
const double bhp_limit,
const double maxPerfPress,
const double vfp_flo_front,
DeferredLogger& deferred_logger) const
{
// Find the bhp-point where production becomes nonzero.
double bhp_max = 0.0;
double low = bhp_limit;
double high = maxPerfPress + 1.0 * unit::barsa;
double f_low = fflo(low);
double f_high = fflo(high);
if constexpr (extraBhpAtThpLimitOutput) {
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + well_.name() +
" low = " + std::to_string(low) +
" high = " + std::to_string(high) +
" f(low) = " + std::to_string(f_low) +
" f(high) = " + std::to_string(f_high));
}
int adjustments = 0;
const int max_adjustments = 10;
const double adjust_amount = 5.0 * unit::barsa;
while (f_low * f_high > 0.0 && adjustments < max_adjustments) {
// Same sign, adjust high to see if we can flip it.
high += adjust_amount;
f_high = fflo(high);
++adjustments;
}
if (f_low * f_high > 0.0) {
if (f_low > 0.0) {
// Even at the BHP limit, we are injecting.
// There will be no solution here, return an
// empty optional.
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE",
"Robust bhp(thp) solve failed due to inoperability for well " + well_.name());
return std::nullopt;
} else {
// Still producing, even at high bhp.
assert(f_high < 0.0);
bhp_max = high;
}
} else {
// Bisect to find a bhp point where we produce, but
// not a large amount ('eps' below).
const double eps = 0.1 * std::fabs(vfp_flo_front);
const int maxit = 50;
int it = 0;
while (std::fabs(f_low) > eps && it < maxit) {
const double curr = 0.5*(low + high);
const double f_curr = fflo(curr);
if (f_curr * f_low > 0.0) {
low = curr;
f_low = f_curr;
} else {
high = curr;
f_high = f_curr;
}
++it;
}
if (it < maxit) {
bhp_max = low;
} else {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE",
"Bisect did not find the bhp-point where we produce for well " + well_.name());
return std::nullopt;
}
}
if constexpr (extraBhpAtThpLimitOutput) {
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + well_.name() +
" low = " + std::to_string(low) +
" high = " + std::to_string(high) +
" f(low) = " + std::to_string(f_low) +
" f(high) = " + std::to_string(f_high) +
" bhp_max = " + std::to_string(bhp_max));
}
return bhp_max;
}
std::optional<double>
WellBhpThpCalculator::
computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& frates,
const std::function<double(const std::vector<double>)>& fbhp,
const std::array<double, 2>& range,
DeferredLogger& deferred_logger) const
{
// Given a VFP function returning bhp as a function of phase
// rates and thp:
// fbhp(rates, thp),
// a function extracting the particular flow rate used for VFP
// lookups:
// flo(rates)
// and the inflow function (assuming the reservoir is fixed):
// frates(bhp)
// we want to solve the equation:
// fbhp(frates(bhp, thplimit)) - bhp = 0
// for bhp.
//
// This may result in 0, 1 or 2 solutions. If two solutions,
// the one corresponding to the lowest bhp (and therefore
// highest rate) should be returned.
// Define the equation we want to solve.
auto eq = [&fbhp, &frates](double bhp) {
return fbhp(frates(bhp)) - bhp;
};
// Find appropriate brackets for the solution.
std::optional<double> approximate_solution;
double low, high;
// trying to use bisect way to locate a bracket
bool finding_bracket = this->bisectBracket(eq, range, low, high, approximate_solution, deferred_logger);
// based on the origional design, if an approximate solution is suggested, we use this value directly
// in the long run, we might change it
if (approximate_solution.has_value()) {
return *approximate_solution;
}
if (!finding_bracket) {
deferred_logger.debug(" Trying the brute force search to bracket the bhp for last attempt ");
finding_bracket = this->bruteForceBracket(eq, range, low, high, deferred_logger);
}
if (!finding_bracket) {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE",
"Robust bhp(thp) solve failed due to not being able to "
"bracket the bhp solution with the brute force search for " + well_.name());
return std::nullopt;
}
// Solve for the proper solution in the given interval.
const int max_iteration = 100;
const double bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0;
try {
const double solved_bhp = RegulaFalsiBisection<ThrowOnError>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
return solved_bhp;
}
catch (...) {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
"Robust bhp(thp) solve failed for well " + well_.name());
return std::nullopt;
}
}
bool
WellBhpThpCalculator::
bisectBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
std::optional<double>& approximate_solution,
DeferredLogger& deferred_logger) const
{
bool finding_bracket = false;
low = range[0];
high = range[1];
double eq_high = eq(high);
double eq_low = eq(low);
const double eq_bhplimit = eq_low;
if constexpr (extraBhpAtThpLimitOutput) {
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + well_.name() +
" low = " + std::to_string(low) +
" high = " + std::to_string(high) +
" eq(low) = " + std::to_string(eq_low) +
" eq(high) = " + std::to_string(eq_high));
}
if (eq_low * eq_high > 0.0) {
// Failed to bracket the zero.
// If this is due to having two solutions, bisect until bracketed.
double abs_low = std::fabs(eq_low);
double abs_high = std::fabs(eq_high);
int bracket_attempts = 0;
const int max_bracket_attempts = 20;
double interval = high - low;
const double min_interval = 1.0 * unit::barsa;
while (eq_low * eq_high > 0.0 && bracket_attempts < max_bracket_attempts && interval > min_interval) {
if (abs_high < abs_low) {
low = 0.5 * (low + high);
eq_low = eq(low);
abs_low = std::fabs(eq_low);
} else {
high = 0.5 * (low + high);
eq_high = eq(high);
abs_high = std::fabs(eq_high);
}
++bracket_attempts;
}
if (eq_low * eq_high <= 0.) {
// We have a bracket!
finding_bracket = true;
// Now, see if (bhplimit, low) is a bracket in addition to (low, high).
// If so, that is the bracket we shall use, choosing the solution with the
// highest flow.
if (eq_low * eq_bhplimit <= 0.0) {
high = low;
low = range[0];
}
} else { // eq_low * eq_high > 0.0
// Still failed bracketing!
const double limit = 0.1 * unit::barsa;
if (std::min(abs_low, abs_high) < limit) {
// 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 " + well_.name());
approximate_solution = abs_low < abs_high ? low : high;
} else {
// Return failure.
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE",
"Robust bhp(thp) solve failed due to bracketing failure for well " +
well_.name());
}
}
} else {
finding_bracket = true;
}
return finding_bracket;
}
bool
WellBhpThpCalculator::
bruteForceBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
DeferredLogger& deferred_logger)
{
bool finding_bracket = false;
low = range[0];
high = range[1];
const int sample_number = 100;
const double interval = (high - low) / sample_number;
double eq_low = eq(low);
double eq_high;
for (int i = 0; i < sample_number + 1; ++i) {
high = range[0] + interval * i;
eq_high = eq(high);
if (eq_high * eq_low <= 0.) {
finding_bracket = true;
break;
}
low = high;
eq_low = eq_high;
}
if (finding_bracket) {
deferred_logger.debug(
" brute force solve found low " + std::to_string(low) + " with eq_low " + std::to_string(eq_low) +
" high " + std::to_string(high) + " with eq_high " + std::to_string(eq_high));
}
return finding_bracket;
}
#define INSTANCE(...) \
template __VA_ARGS__ WellBhpThpCalculator:: \
calculateBhpFromThp<__VA_ARGS__>(const WellState&, \
const std::vector<__VA_ARGS__>&, \
const Well&, \
const SummaryState&, \
const double, \
DeferredLogger&) const;
INSTANCE(double)
INSTANCE(DenseAd::Evaluation<double,3,0u>)
INSTANCE(DenseAd::Evaluation<double,4,0u>)
INSTANCE(DenseAd::Evaluation<double,5,0u>)
INSTANCE(DenseAd::Evaluation<double,6,0u>)
INSTANCE(DenseAd::Evaluation<double,7,0u>)
INSTANCE(DenseAd::Evaluation<double,8,0u>)
INSTANCE(DenseAd::Evaluation<double,9,0u>)
INSTANCE(DenseAd::Evaluation<double,-1,4u>)
INSTANCE(DenseAd::Evaluation<double,-1,5u>)
INSTANCE(DenseAd::Evaluation<double,-1,6u>)
INSTANCE(DenseAd::Evaluation<double,-1,7u>)
INSTANCE(DenseAd::Evaluation<double,-1,8u>)
INSTANCE(DenseAd::Evaluation<double,-1,9u>)
INSTANCE(DenseAd::Evaluation<double,-1,10u>)
} // namespace Opm
opm/simulators/wells/WellBhpThpCalculator.hpp
+141
View File
@@ -0,0 +1,141 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2017 IRIS
Copyright 2019 Norce
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELL_BPH_THP_CALCULATOR_HEADER_INCLUDED
#define OPM_WELL_BPH_THP_CALCULATOR_HEADER_INCLUDED
#include <functional>
#include <optional>
#include <string>
#include <vector>
namespace Opm
{
class DeferredLogger;
class SummaryState;
class Well;
class WellInterfaceGeneric;
class WellState;
//! \brief Class for computing BHP limits.
class WellBhpThpCalculator {
public:
//! \brief Constructor sets reference to well.
WellBhpThpCalculator(const WellInterfaceGeneric& well) : well_(well) {}
//! \brief Checks if well has THP constraints.
bool wellHasTHPConstraints(const SummaryState& summaryState) const;
//! \brief Get THP constraint for well.
double getTHPConstraint(const SummaryState& summaryState) const;
//! \brief Obtain the most strict BHP from BHP limits.
double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
//! \brief Calculates THP from BHP.
double calculateThpFromBhp(const std::vector<double>& rates,
const double bhp,
const double rho,
const double alq,
DeferredLogger& deferred_logger) const;
//! \brief Compute BHP from THP limit for a producer.
std::optional<double>
computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double maxPerfPress,
const double rho,
const double alq_value,
DeferredLogger& deferred_logger) const;
//! \brief Compute BHP from THP limit for an injector.
std::optional<double>
computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double rho,
const double flo_rel_tol,
const int max_iteration,
const bool throwOnError,
DeferredLogger& deferred_logger) const;
//! \brief Update THP.
void updateThp(const double rho,
const std::function<double()>& alq_value,
const std::array<unsigned,3>& active,
WellState& well_state,
DeferredLogger& deferred_logger) const;
template<class EvalWell>
EvalWell calculateBhpFromThp(const WellState& well_state,
const std::vector<EvalWell>& rates,
const Well& well,
const SummaryState& summaryState,
const double rho,
DeferredLogger& deferred_logger) const;
private:
//! \brief Compute BHP from THP limit for an injector - implementation.
template<class ErrorPolicy>
std::optional<double>
computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double rho,
const double flo_rel_tol,
const int max_iteration,
DeferredLogger& deferred_logger) const;
//! \brief Calculate max BHP.
std::optional<double>
bhpMax(const std::function<double(const double)>& fflo,
const double bhp_limit,
const double maxPerfPress,
const double vfp_flo_front,
DeferredLogger& deferred_logger) const;
//! \brief Common code for finding BHP from THP limit for producers/injectors.
std::optional<double>
computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& frates,
const std::function<double(const std::vector<double>)>& fbhp,
const std::array<double, 2>& range,
DeferredLogger& deferred_logger) const;
//! \brief Find limits using bisection.
bool bisectBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
std::optional<double>& approximate_solution,
DeferredLogger& deferred_logger) const;
//! \brief Find limits using brute-force solver.
static bool bruteForceBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
DeferredLogger& deferred_logger);
const WellInterfaceGeneric& well_; //!< Reference to well interface
};
}
#endif // OPM_WELL_BHP_THP_CALCULATOR_HEADER_INCLUDED
opm/simulators/wells/WellInterface.hpp
+2 -2
View File
@@ -171,8 +171,8 @@ public:
virtual std::optional<double> computeBhpAtThpLimitProdWithAlq(
const Simulator& ebos_simulator,
const SummaryState& summary_state,
DeferredLogger& deferred_logger,
double alq_value
const double alq_value,
DeferredLogger& deferred_logger
) const = 0;
/// using the solution x to recover the solution xw for wells and applying
opm/simulators/wells/WellInterfaceEval.cpp
+1 -65
View File
@@ -260,52 +260,6 @@ assembleControlEqInj_(const WellState& well_state,
}
}
template<class FluidSystem>
template<class EvalWell>
EvalWell
WellInterfaceEval<FluidSystem>::
calculateBhpFromThp(const WellState& well_state,
const std::vector<EvalWell>& rates,
const Well& well,
const SummaryState& summaryState,
const double rho,
DeferredLogger& deferred_logger) const
{
// TODO: when well is under THP control, the BHP is dependent on the rates,
// the well rates is also dependent on the BHP, so it might need to do some iteration.
// However, when group control is involved, change of the rates might impacts other wells
// so iterations on a higher level will be required. Some investigation might be needed when
// we face problems under THP control.
assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
const EvalWell aqua = rates[Water];
const EvalWell liquid = rates[Oil];
const EvalWell vapour = rates[Gas];
// pick the reference density
// typically the reference in the top layer
if (baseif_.isInjector() )
{
const auto& controls = well.injectionControls(summaryState);
const double vfp_ref_depth = baseif_.vfpProperties()->getInj()->getTable(controls.vfp_table_number).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
return baseif_.vfpProperties()->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, baseif_.getTHPConstraint(summaryState)) - dp;
}
else if (baseif_.isProducer()) {
const auto& controls = well.productionControls(summaryState);
const double vfp_ref_depth = baseif_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
const auto& wfr = baseif_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, baseif_.indexOfWell());
const auto& gfr = baseif_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, baseif_.indexOfWell());
const bool use_vfpexplicit = baseif_.useVfpExplicit();
return baseif_.vfpProperties()->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, baseif_.getTHPConstraint(summaryState), baseif_.getALQ(well_state), wfr, gfr, use_vfpexplicit) - dp;
}
else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + baseif_.name(), deferred_logger);
}
}
#define INSTANCE_METHODS(A,...) \
template void WellInterfaceEval<A>:: \
assembleControlEqProd_<__VA_ARGS__>(const WellState&, \
@@ -328,14 +282,7 @@ assembleControlEqInj_<__VA_ARGS__>(const WellState&, \
const __VA_ARGS__&, \
const std::function<__VA_ARGS__()>&, \
__VA_ARGS__&, \
DeferredLogger&) const; \
template __VA_ARGS__ WellInterfaceEval<A>:: \
calculateBhpFromThp<__VA_ARGS__>(const WellState&, \
const std::vector<__VA_ARGS__>&, \
const Well&, \
const SummaryState&, \
const double, \
DeferredLogger&) const;
DeferredLogger&) const;
using FluidSys = BlackOilFluidSystem<double, BlackOilDefaultIndexTraits>;
@@ -356,15 +303,4 @@ INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,8u>)
INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,9u>)
INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,10u>)
#define INSTANCE_BHP(...) \
template double WellInterfaceEval<__VA_ARGS__>:: \
calculateBhpFromThp<double>(const WellState&, \
const std::vector<double>&, \
const Well&, \
const SummaryState&, \
const double, \
DeferredLogger&) const;
INSTANCE_BHP(FluidSys)
} // namespace Opm
opm/simulators/wells/WellInterfaceEval.hpp
-8
View File
@@ -49,14 +49,6 @@ class WellInterfaceEval {
static constexpr int Gas = BlackoilPhases::Vapour;
public:
template <class EvalWell>
EvalWell calculateBhpFromThp(const WellState& well_state,
const std::vector<EvalWell>& rates,
const Well& well,
const SummaryState& summaryState,
const double rho,
DeferredLogger& deferred_logger) const;
template<class EvalWell, class BhpFromThpFunc>
void assembleControlEqProd(const WellState& well_state,
const GroupState& group_state,
opm/simulators/wells/WellInterfaceGeneric.cpp
+3 -361
View File
@@ -23,12 +23,12 @@
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
#include <opm/common/utility/numeric/RootFinders.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/PerforationData.hpp>
#include <opm/simulators/wells/ParallelWellInfo.hpp>
#include <opm/simulators/wells/VFPHelpers.hpp>
#include <opm/simulators/wells/VFPProperties.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/wells/WellHelpers.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <opm/simulators/wells/WellTest.hpp>
@@ -185,35 +185,7 @@ bool WellInterfaceGeneric::wellHasTHPConstraints(const SummaryState& summaryStat
return true;
}
if (well_ecl_.isInjector()) {
const auto controls = well_ecl_.injectionControls(summaryState);
if (controls.hasControl(Well::InjectorCMode::THP))
return true;
}
if (well_ecl_.isProducer( )) {
const auto controls = well_ecl_.productionControls(summaryState);
if (controls.hasControl(Well::ProducerCMode::THP))
return true;
}
return false;
}
double WellInterfaceGeneric::mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const
{
if (well_ecl_.isInjector()) {
const auto& controls = well_ecl_.injectionControls(summaryState);
return controls.bhp_limit;
}
if (well_ecl_.isProducer( )) {
const auto& controls = well_ecl_.productionControls(summaryState);
return controls.bhp_limit;
}
return 0.0;
return WellBhpThpCalculator(*this).wellHasTHPConstraints(summaryState);
}
void WellInterfaceGeneric::updateWellTestState(const SingleWellState& ws,
@@ -233,23 +205,13 @@ void WellInterfaceGeneric::updateWellTestState(const SingleWellState& ws,
// TODO: well can be shut/closed due to other reasons
}
double WellInterfaceGeneric::getTHPConstraint(const SummaryState& summaryState) const
{
if (dynamic_thp_limit_) {
return *dynamic_thp_limit_;
}
if (well_ecl_.isInjector()) {
const auto& controls = well_ecl_.injectionControls(summaryState);
return controls.thp_limit;
}
if (well_ecl_.isProducer( )) {
const auto& controls = well_ecl_.productionControls(summaryState);
return controls.thp_limit;
}
return 0.0;
return WellBhpThpCalculator(*this).getTHPConstraint(summaryState);
}
bool WellInterfaceGeneric::underPredictionMode() const
@@ -454,324 +416,4 @@ void WellInterfaceGeneric::reportWellSwitching(const SingleWellState& ws, Deferr
}
}
std::optional<double>
WellInterfaceGeneric::
bhpMax(const std::function<double(const double)>& fflo,
const double bhp_limit,
const double maxPerfPress,
const double vfp_flo_front,
DeferredLogger& deferred_logger) const
{
// Find the bhp-point where production becomes nonzero.
double bhp_max = 0.0;
double low = bhp_limit;
double high = maxPerfPress + 1.0 * unit::barsa;
double f_low = fflo(low);
double f_high = fflo(high);
if constexpr (extraBhpAtThpLimitProdOutput) {
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + this->name() +
" low = " + std::to_string(low) +
" high = " + std::to_string(high) +
" f(low) = " + std::to_string(f_low) +
" f(high) = " + std::to_string(f_high));
}
int adjustments = 0;
const int max_adjustments = 10;
const double adjust_amount = 5.0 * unit::barsa;
while (f_low * f_high > 0.0 && adjustments < max_adjustments) {
// Same sign, adjust high to see if we can flip it.
high += adjust_amount;
f_high = fflo(high);
++adjustments;
}
if (f_low * f_high > 0.0) {
if (f_low > 0.0) {
// Even at the BHP limit, we are injecting.
// There will be no solution here, return an
// empty optional.
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE",
"Robust bhp(thp) solve failed due to inoperability for well " + this->name());
return std::nullopt;
} else {
// Still producing, even at high bhp.
assert(f_high < 0.0);
bhp_max = high;
}
} else {
// Bisect to find a bhp point where we produce, but
// not a large amount ('eps' below).
const double eps = 0.1 * std::fabs(vfp_flo_front);
const int maxit = 50;
int it = 0;
while (std::fabs(f_low) > eps && it < maxit) {
const double curr = 0.5*(low + high);
const double f_curr = fflo(curr);
if (f_curr * f_low > 0.0) {
low = curr;
f_low = f_curr;
} else {
high = curr;
f_high = f_curr;
}
++it;
}
if (it < maxit) {
bhp_max = low;
} else {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE",
"Bisect did not find the bhp-point where we produce for well " + this->name());
return std::nullopt;
}
}
if constexpr (extraBhpAtThpLimitProdOutput) {
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + this->name() +
" low = " + std::to_string(low) +
" high = " + std::to_string(high) +
" f(low) = " + std::to_string(f_low) +
" f(high) = " + std::to_string(f_high) +
" bhp_max = " + std::to_string(bhp_max));
}
return bhp_max;
}
bool
WellInterfaceGeneric::
bisectBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
std::optional<double>& approximate_solution,
DeferredLogger& deferred_logger) const
{
bool finding_bracket = false;
low = range[0];
high = range[1];
double eq_high = eq(high);
double eq_low = eq(low);
const double eq_bhplimit = eq_low;
if constexpr (extraBhpAtThpLimitProdOutput) {
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + this->name() +
" low = " + std::to_string(low) +
" high = " + std::to_string(high) +
" eq(low) = " + std::to_string(eq_low) +
" eq(high) = " + std::to_string(eq_high));
}
if (eq_low * eq_high > 0.0) {
// Failed to bracket the zero.
// If this is due to having two solutions, bisect until bracketed.
double abs_low = std::fabs(eq_low);
double abs_high = std::fabs(eq_high);
int bracket_attempts = 0;
const int max_bracket_attempts = 20;
double interval = high - low;
const double min_interval = 1.0 * unit::barsa;
while (eq_low * eq_high > 0.0 && bracket_attempts < max_bracket_attempts && interval > min_interval) {
if (abs_high < abs_low) {
low = 0.5 * (low + high);
eq_low = eq(low);
abs_low = std::fabs(eq_low);
} else {
high = 0.5 * (low + high);
eq_high = eq(high);
abs_high = std::fabs(eq_high);
}
++bracket_attempts;
}
if (eq_low * eq_high <= 0.) {
// We have a bracket!
finding_bracket = true;
// Now, see if (bhplimit, low) is a bracket in addition to (low, high).
// If so, that is the bracket we shall use, choosing the solution with the
// highest flow.
if (eq_low * eq_bhplimit <= 0.0) {
high = low;
low = range[0];
}
} else { // eq_low * eq_high > 0.0
// Still failed bracketing!
const double limit = 0.1 * unit::barsa;
if (std::min(abs_low, abs_high) < limit) {
// 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;
} else {
// Return failure.
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE",
"Robust bhp(thp) solve failed due to bracketing failure for well " +
this->name());
}
}
} else {
finding_bracket = true;
}
return finding_bracket;
}
bool
WellInterfaceGeneric::
bruteForceBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
DeferredLogger& deferred_logger) const
{
bool finding_bracket = false;
low = range[0];
high = range[1];
const int sample_number = 100;
const double interval = (high - low) / sample_number;
double eq_low = eq(low);
double eq_high;
for (int i = 0; i < sample_number + 1; ++i) {
high = range[0] + interval * i;
eq_high = eq(high);
if (eq_high * eq_low <= 0.) {
finding_bracket = true;
break;
}
low = high;
eq_low = eq_high;
}
if (finding_bracket) {
deferred_logger.debug(
" brute force solve found low " + std::to_string(low) + " with eq_low " + std::to_string(eq_low) +
" high " + std::to_string(high) + " with eq_high " + std::to_string(eq_high));
}
return finding_bracket;
}
std::optional<double>
WellInterfaceGeneric::
computeBhpAtThpLimitProdCommon(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double maxPerfPress,
const double rho,
const double alq_value,
DeferredLogger& deferred_logger) const
{
// Given a VFP function returning bhp as a function of phase
// rates and thp:
// fbhp(rates, thp),
// a function extracting the particular flow rate used for VFP
// lookups:
// flo(rates)
// and the inflow function (assuming the reservoir is fixed):
// frates(bhp)
// we want to solve the equation:
// fbhp(frates(bhp, thplimit)) - bhp = 0
// for bhp.
//
// This may result in 0, 1 or 2 solutions. If two solutions,
// the one corresponding to the lowest bhp (and therefore
// highest rate) should be returned.
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
// Make the fbhp() function.
const auto& controls = this->wellEcl().productionControls(summary_state);
const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const double vfp_ref_depth = table.getDatumDepth();
const double thp_limit = this->getTHPConstraint(summary_state);
const double dp = wellhelpers::computeHydrostaticCorrection(this->refDepth(), vfp_ref_depth, rho, this->gravity());
auto fbhp = [this, &controls, thp_limit, dp, alq_value](const std::vector<double>& rates) {
assert(rates.size() == 3);
const auto& wfr = this->vfpProperties()->getExplicitWFR(controls.vfp_table_number, this->indexOfWell());
const auto& gfr = this->vfpProperties()->getExplicitGFR(controls.vfp_table_number, this->indexOfWell());
const bool use_vfpexp = this->useVfpExplicit();
return this->vfpProperties()->getProd()
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, alq_value, wfr, gfr, use_vfpexp) - dp;
};
// Make the flo() function.
auto flo = [&table](const std::vector<double>& rates) {
return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]);
};
// Find the bhp-point where production becomes nonzero.
auto fflo = [&flo, &frates](double bhp) { return flo(frates(bhp)); };
auto bhp_max = this->bhpMax(fflo, controls.bhp_limit, maxPerfPress, table.getFloAxis().front(), deferred_logger);
// could not solve for the bhp-point, we could not continue to find the bhp
if (!bhp_max.has_value()) {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE",
"Robust bhp(thp) solve failed due to not being able to "
"find bhp-point where production becomes non-zero for well " + this->name());
return std::nullopt;
}
const std::array<double, 2> range {controls.bhp_limit, *bhp_max};
return computeBhpAtThpLimitCommon(frates, fbhp, range, deferred_logger);
}
std::optional<double>
WellInterfaceGeneric::
computeBhpAtThpLimitCommon(const std::function<std::vector<double>(const double)>& frates,
const std::function<double(const std::vector<double>)>& fbhp,
const std::array<double, 2>& range,
DeferredLogger& deferred_logger) const
{
// Given a VFP function returning bhp as a function of phase
// rates and thp:
// fbhp(rates, thp),
// a function extracting the particular flow rate used for VFP
// lookups:
// flo(rates)
// and the inflow function (assuming the reservoir is fixed):
// frates(bhp)
// we want to solve the equation:
// fbhp(frates(bhp, thplimit)) - bhp = 0
// for bhp.
//
// This may result in 0, 1 or 2 solutions. If two solutions,
// the one corresponding to the lowest bhp (and therefore
// highest rate) should be returned.
// Define the equation we want to solve.
auto eq = [&fbhp, &frates](double bhp) {
return fbhp(frates(bhp)) - bhp;
};
// Find appropriate brackets for the solution.
std::optional<double> approximate_solution;
double low, high;
// trying to use bisect way to locate a bracket
bool finding_bracket = this->bisectBracket(eq, range, low, high, approximate_solution, deferred_logger);
// based on the origional design, if an approximate solution is suggested, we use this value directly
// in the long run, we might change it
if (approximate_solution.has_value()) {
return *approximate_solution;
}
if (!finding_bracket) {
deferred_logger.debug(" Trying the brute force search to bracket the bhp for last attempt ");
finding_bracket = this->bruteForceBracket(eq, range, low, high, deferred_logger);
}
if (!finding_bracket) {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE",
"Robust bhp(thp) solve failed due to not being able to "
"bracket the bhp solution with the brute force search for " + this->name());
return std::nullopt;
}
// Solve for the proper solution in the given interval.
const int max_iteration = 100;
const double bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0;
try {
const double solved_bhp = RegulaFalsiBisection<ThrowOnError>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
return solved_bhp;
}
catch (...) {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
"Robust bhp(thp) solve failed for well " + this->name());
return std::nullopt;
}
}
} // namespace Opm
opm/simulators/wells/WellInterfaceGeneric.hpp
-36
View File
@@ -50,9 +50,6 @@ class Schedule;
class WellInterfaceGeneric {
public:
static constexpr bool extraBhpAtThpLimitProdOutput = false;
WellInterfaceGeneric(const Well& well,
const ParallelWellInfo& parallel_well_info,
const int time_step,
@@ -189,13 +186,6 @@ public:
bool changedToOpenThisStep() const {
return this->changed_to_open_this_step_;
}
std::optional<double> computeBhpAtThpLimitProdCommon(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double maxPerfPress,
const double rho,
const double alq_value,
DeferredLogger& deferred_logger
) const;
void updateWellTestState(const SingleWellState& ws,
const double& simulationTime,
@@ -205,32 +195,6 @@ public:
protected:
bool getAllowCrossFlow() const;
double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
std::optional<double> bhpMax(const std::function<double(const double)>& fflo,
const double bhp_limit,
const double maxPerfPress,
const double vfp_flo_front,
DeferredLogger& deferred_logger) const;
std::optional<double> computeBhpAtThpLimitCommon(
const std::function<std::vector<double>(const double)>& frates,
const std::function<double(const std::vector<double>)>& fbhp,
const std::array<double, 2>& range,
DeferredLogger& deferred_logger) const;
bool bruteForceBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
DeferredLogger& deferred_logger) const;
bool bisectBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
std::optional<double>& approximate_solution,
DeferredLogger& deferred_logger) const;
// definition of the struct OperabilityStatus
struct OperabilityStatus {
opm/simulators/wells/WellInterface_impl.hpp
+13 -2
View File
@@ -23,6 +23,7 @@
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/GroupState.hpp>
#include <opm/simulators/wells/TargetCalculator.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <dune/common/version.hh>
@@ -811,7 +812,12 @@ namespace Opm
case Well::InjectorCMode::THP:
{
auto rates = ws.surface_rates;
double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
rates,
well,
summaryState,
this->getRefDensity(),
deferred_logger);
ws.bhp = bhp;
ws.thp = this->getTHPConstraint(summaryState);
@@ -1036,7 +1042,12 @@ namespace Opm
{
auto rates = ws.surface_rates;
this->adaptRatesForVFP(rates);
double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
double bhp = WellBhpThpCalculator(*this).calculateBhpFromThp(well_state,
rates,
well,
summaryState,
this->getRefDensity(),
deferred_logger);
ws.bhp = bhp;
ws.thp = this->getTHPConstraint(summaryState);