OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

changed: move computeBhpFromThpLimitInj to WellBhpThpCalculator

Browse Source
This commit is contained in:
Arne Morten Kvarving 2022-10-19 09:55:14 +02:00
parent 6e214557e2
commit f214ccc138
9 changed files with 253 additions and 394 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

189
opm/simulators/wells/MultisegmentWellGeneric.cpp
View File

@ -207,195 +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>::

5
opm/simulators/wells/MultisegmentWellGeneric.hpp
View File

@ -60,11 +60,6 @@ protected:
/// number of segments for this well
int numberOfSegments() 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,

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

@ -1953,13 +1953,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) {
@ -1971,8 +1974,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);
}

188
opm/simulators/wells/StandardWellGeneric.cpp
View File

@ -167,194 +167,6 @@ computeBhpAtThpLimitProdWithAlq(const std::function<std::vector<double>(const do
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>::

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

@ -81,9 +81,6 @@ protected:
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,

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

@ -2655,9 +2655,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);
}

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

@ -36,6 +36,7 @@
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <cassert>
#include <cmath>
static constexpr bool extraBhpAtThpLimitOutput = false;
@ -197,6 +198,217 @@ computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>&
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);
}
}
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,

20
opm/simulators/wells/WellBhpThpCalculator.hpp
View File

@ -67,7 +67,27 @@ public:
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;
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,

1
opm/simulators/wells/WellInterfaceGeneric.cpp
View File

@ -23,7 +23,6 @@
#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>