mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Move computeBhpAtThpLimitProd MSW code to common and apply it for both MSW and STW
The computeBhpAtThpLimitProd for MSW is faster and sufficiently robust and should also be used by STW
This commit is contained in:
Tor Harald Sandve
@@ -444,284 +444,7 @@ computeBhpAtThpLimitProdWithAlq(
|
|||||||
DeferredLogger& deferred_logger,
|
DeferredLogger& deferred_logger,
|
||||||
double alq_value) const
|
double alq_value) const
|
||||||
{
|
{
|
||||||
// Given a VFP function returning bhp as a function of phase
|
return baseif_.computeBhpAtThpLimitProdCommon(frates, summary_state, maxPerfPress, rho, alq_value, deferred_logger);
|
||||||
// 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 = baseif_.wellEcl().productionControls(summary_state);
|
|
||||||
const auto& table = baseif_.vfpProperties()->getProd()->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, alq_value](
|
|
||||||
const std::vector<double>& rates) {
|
|
||||||
assert(rates.size() == 3);
|
|
||||||
return baseif_.vfpProperties()->getProd()
|
|
||||||
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, alq_value) - 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 " + baseif_.name());
|
|
||||||
return std::nullopt;
|
|
||||||
}
|
|
||||||
// Define the equation we want to solve.
|
|
||||||
auto eq = [&fbhp, &frates](double bhp) {
|
|
||||||
return fbhp(frates(bhp)) - bhp;
|
|
||||||
};
|
|
||||||
|
|
||||||
// Find appropriate brackets for the solution.
|
|
||||||
const std::array<double, 2> range {controls.bhp_limit, *bhp_max};
|
|
||||||
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 " + baseif_.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 " + baseif_.name());
|
|
||||||
return std::nullopt;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
template<typename Scalar>
|
|
||||||
std::optional<double>
|
|
||||||
MultisegmentWellGeneric<Scalar>::
|
|
||||||
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);
|
|
||||||
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + baseif_.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 " + baseif_.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 " + baseif_.name());
|
|
||||||
return std::nullopt;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + baseif_.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;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
template<typename Scalar>
|
|
||||||
bool
|
|
||||||
MultisegmentWellGeneric<Scalar>::
|
|
||||||
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;
|
|
||||||
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + baseif_.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 = 3.0 * unit::barsa;
|
|
||||||
if (std::min(abs_low, abs_high) < limit) {
|
|
||||||
// Return the least bad solution if less off than 3 bar.
|
|
||||||
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE",
|
|
||||||
"Robust bhp(thp) not solved precisely for well " + baseif_.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 " +
|
|
||||||
baseif_.name());
|
|
||||||
}
|
|
||||||
}
|
|
||||||
} else {
|
|
||||||
finding_bracket = true;
|
|
||||||
}
|
|
||||||
return finding_bracket;
|
|
||||||
}
|
|
||||||
|
|
||||||
template<typename Scalar>
|
|
||||||
bool
|
|
||||||
MultisegmentWellGeneric<Scalar>::
|
|
||||||
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;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
template<typename Scalar>
|
template<typename Scalar>
|
||||||
|
|||||||
@@ -199,8 +199,13 @@ StandardWellGeneric<Scalar>::
|
|||||||
computeBhpAtThpLimitProdWithAlq(const std::function<std::vector<double>(const double)>& frates,
|
computeBhpAtThpLimitProdWithAlq(const std::function<std::vector<double>(const double)>& frates,
|
||||||
const SummaryState& summary_state,
|
const SummaryState& summary_state,
|
||||||
DeferredLogger& deferred_logger,
|
DeferredLogger& deferred_logger,
|
||||||
|
double maxPerfPress,
|
||||||
double alq_value) const
|
double alq_value) const
|
||||||
{
|
{
|
||||||
|
|
||||||
|
if (true) {
|
||||||
|
return baseif_.computeBhpAtThpLimitProdCommon(frates, summary_state, maxPerfPress, getRho(), alq_value, deferred_logger);
|
||||||
|
}
|
||||||
// Given a VFP function returning bhp as a function of phase
|
// Given a VFP function returning bhp as a function of phase
|
||||||
// rates and thp:
|
// rates and thp:
|
||||||
// fbhp(rates, thp),
|
// fbhp(rates, thp),
|
||||||
|
|||||||
@@ -105,6 +105,7 @@ protected:
|
|||||||
std::optional<double> computeBhpAtThpLimitProdWithAlq(const std::function<std::vector<double>(const double)>& frates,
|
std::optional<double> computeBhpAtThpLimitProdWithAlq(const std::function<std::vector<double>(const double)>& frates,
|
||||||
const SummaryState& summary_state,
|
const SummaryState& summary_state,
|
||||||
DeferredLogger& deferred_logger,
|
DeferredLogger& deferred_logger,
|
||||||
|
double maxPerfPress,
|
||||||
double alq_value) const;
|
double alq_value) const;
|
||||||
|
|
||||||
// Base interface reference
|
// Base interface reference
|
||||||
|
|||||||
@@ -2320,9 +2320,18 @@ namespace Opm
|
|||||||
return rates;
|
return rates;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
double max_pressure = 0.0;
|
||||||
|
for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
|
||||||
|
const int cell_idx = this->well_cells_[perf];
|
||||||
|
const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
|
||||||
|
const auto& fs = int_quants.fluidState();
|
||||||
|
double pressure_cell = this->getPerfCellPressure(fs).value();
|
||||||
|
max_pressure = std::max(max_pressure, pressure_cell);
|
||||||
|
}
|
||||||
return this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitProdWithAlq(frates,
|
return this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitProdWithAlq(frates,
|
||||||
summary_state,
|
summary_state,
|
||||||
deferred_logger,
|
deferred_logger,
|
||||||
|
max_pressure,
|
||||||
alq_value);
|
alq_value);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|||||||
@@ -23,12 +23,14 @@
|
|||||||
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
|
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
|
||||||
|
|
||||||
#include <opm/input/eclipse/Schedule/Well/WellTestState.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/utils/DeferredLoggingErrorHelpers.hpp>
|
||||||
#include <opm/simulators/wells/PerforationData.hpp>
|
#include <opm/simulators/wells/PerforationData.hpp>
|
||||||
#include <opm/simulators/wells/ParallelWellInfo.hpp>
|
#include <opm/simulators/wells/ParallelWellInfo.hpp>
|
||||||
#include <opm/simulators/wells/VFPProperties.hpp>
|
#include <opm/simulators/wells/VFPProperties.hpp>
|
||||||
#include <opm/simulators/wells/WellState.hpp>
|
#include <opm/simulators/wells/WellState.hpp>
|
||||||
|
#include <opm/simulators/wells/WellHelpers.hpp>
|
||||||
|
#include <opm/simulators/wells/VFPHelpers.hpp>
|
||||||
#include <cassert>
|
#include <cassert>
|
||||||
#include <cmath>
|
#include <cmath>
|
||||||
#include <cstddef>
|
#include <cstddef>
|
||||||
@@ -408,5 +410,288 @@ void WellInterfaceGeneric::reportWellSwitching(const SingleWellState& ws, Deferr
|
|||||||
deferred_logger.info(msg);
|
deferred_logger.info(msg);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
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);
|
||||||
|
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;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
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;
|
||||||
|
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 = 3.0 * unit::barsa;
|
||||||
|
if (std::min(abs_low, abs_high) < limit) {
|
||||||
|
// Return the least bad solution if less off than 3 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);
|
||||||
|
return this->vfpProperties()->getProd()
|
||||||
|
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, alq_value) - 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;
|
||||||
|
}
|
||||||
|
// Define the equation we want to solve.
|
||||||
|
auto eq = [&fbhp, &frates](double bhp) {
|
||||||
|
return fbhp(frates(bhp)) - bhp;
|
||||||
|
};
|
||||||
|
|
||||||
|
// Find appropriate brackets for the solution.
|
||||||
|
const std::array<double, 2> range {controls.bhp_limit, *bhp_max};
|
||||||
|
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
|
} // namespace Opm
|
||||||
|
|||||||
@@ -177,6 +177,13 @@ public:
|
|||||||
bool changedToOpenThisStep() const {
|
bool changedToOpenThisStep() const {
|
||||||
return this->changed_to_open_this_step_;
|
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;
|
||||||
protected:
|
protected:
|
||||||
bool getAllowCrossFlow() const;
|
bool getAllowCrossFlow() const;
|
||||||
double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
|
double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
|
||||||
@@ -185,6 +192,23 @@ protected:
|
|||||||
WellTestState& well_test_state,
|
WellTestState& well_test_state,
|
||||||
DeferredLogger& deferred_logger) const;
|
DeferredLogger& deferred_logger) 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;
|
||||||
|
|
||||||
|
|
||||||
// definition of the struct OperabilityStatus
|
// definition of the struct OperabilityStatus
|
||||||
struct OperabilityStatus {
|
struct OperabilityStatus {
|
||||||
|
|||||||
Reference in New Issue
Block a user