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

WellBhpThpCalculator: template Scalar type

Browse Source
This commit is contained in:
Arne Morten Kvarving 2024-02-20 11:14:54 +01:00
parent 5a340258aa
commit 9c431d1921
2 changed files with 296 additions and 264 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -48,11 +48,11 @@
static constexpr bool extraBhpAtThpLimitOutput = false; static constexpr bool extraBhpAtThpLimitOutput = false;
static constexpr bool extraThpFromBhpOutput = false; static constexpr bool extraThpFromBhpOutput = false;
namespace Opm namespace Opm {
{
bool template<class Scalar>
WellBhpThpCalculator::wellHasTHPConstraints(const SummaryState& summaryState) const bool WellBhpThpCalculator<Scalar>::
wellHasTHPConstraints(const SummaryState& summaryState) const
{ {
const auto& well_ecl = well_.wellEcl(); const auto& well_ecl = well_.wellEcl();
if (well_ecl.isInjector()) { if (well_ecl.isInjector()) {
@ -70,7 +70,9 @@ WellBhpThpCalculator::wellHasTHPConstraints(const SummaryState& summaryState) co
return false; return false;
} }
double WellBhpThpCalculator::getTHPConstraint(const SummaryState& summaryState) const template<class Scalar>
Scalar WellBhpThpCalculator<Scalar>::
getTHPConstraint(const SummaryState& summaryState) const
{ {
const auto& well_ecl = well_.wellEcl(); const auto& well_ecl = well_.wellEcl();
if (well_ecl.isInjector()) { if (well_ecl.isInjector()) {
@ -86,7 +88,9 @@ double WellBhpThpCalculator::getTHPConstraint(const SummaryState& summaryState)
return 0.0; return 0.0;
} }
double WellBhpThpCalculator::mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const template<class Scalar>
Scalar WellBhpThpCalculator<Scalar>::
mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const
{ {
const auto& well_ecl = well_.wellEcl(); const auto& well_ecl = well_.wellEcl();
if (well_ecl.isInjector()) { if (well_ecl.isInjector()) {
@ -102,12 +106,14 @@ double WellBhpThpCalculator::mostStrictBhpFromBhpLimits(const SummaryState& summ
return 0.0; return 0.0;
} }
double WellBhpThpCalculator::calculateThpFromBhp(const std::vector<double>& rates, template<class Scalar>
const double bhp, Scalar WellBhpThpCalculator<Scalar>::
const double rho, calculateThpFromBhp(const std::vector<Scalar>& rates,
const std::optional<double>& alq, const Scalar bhp,
const double thp_limit, const Scalar rho,
DeferredLogger& deferred_logger) const const std::optional<Scalar>& alq,
const Scalar thp_limit,
DeferredLogger& deferred_logger) const
{ {
assert(int(rates.size()) == 3); // the vfp related only supports three phases now. assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
@ -115,31 +121,31 @@ double WellBhpThpCalculator::calculateThpFromBhp(const std::vector<double>& rate
static constexpr int Oil = BlackoilPhases::Liquid; static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour; static constexpr int Gas = BlackoilPhases::Vapour;
const double aqua = rates[Water]; const Scalar aqua = rates[Water];
const double liquid = rates[Oil]; const Scalar liquid = rates[Oil];
const double vapour = rates[Gas]; const Scalar vapour = rates[Gas];
// pick the density in the top layer // pick the density in the top layer
double thp = 0.0; Scalar thp = 0.0;
const int table_id = well_.wellEcl().vfp_table_number(); const int table_id = well_.wellEcl().vfp_table_number();
if (well_.isInjector()) { if (well_.isInjector()) {
assert(!alq.has_value()); assert(!alq.has_value());
const double vfp_ref_depth = well_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth(); const Scalar vfp_ref_depth = well_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
auto thp_func = auto thp_func =
[this, table_id, aqua, liquid, vapour, dp]( [this, table_id, aqua, liquid, vapour, dp](
const double bhp_value, const double pressure_loss) { const Scalar bhp_value, const Scalar pressure_loss) {
return this->well_.vfpProperties()->getInj()->thp( return this->well_.vfpProperties()->getInj()->thp(
table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss); table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss);
}; };
thp = findThpFromBhpIteratively(thp_func, bhp, thp_limit, dp, deferred_logger); thp = findThpFromBhpIteratively(thp_func, bhp, thp_limit, dp, deferred_logger);
} }
else if (well_.isProducer()) { else if (well_.isProducer()) {
const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth(); const Scalar vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
auto thp_func = auto thp_func =
[this, table_id, aqua, liquid, vapour, dp, &alq] [this, table_id, aqua, liquid, vapour, dp, &alq]
(const double bhp_value, const double pressure_loss) { (const Scalar bhp_value, const Scalar pressure_loss) {
return this->well_.vfpProperties()->getProd()->thp( return this->well_.vfpProperties()->getProd()->thp(
table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss, alq.value()); table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss, alq.value());
}; };
@ -151,36 +157,35 @@ double WellBhpThpCalculator::calculateThpFromBhp(const std::vector<double>& rate
return thp; return thp;
} }
double template<class Scalar>
WellBhpThpCalculator:: Scalar WellBhpThpCalculator<Scalar>::
findThpFromBhpIteratively( findThpFromBhpIteratively(const std::function<Scalar(const Scalar, const Scalar)>& thp_func,
const std::function<double(const double, const double)>& thp_func, const Scalar bhp,
const double bhp, const Scalar thp_limit,
const double thp_limit, const Scalar dp,
const double dp, DeferredLogger& deferred_logger) const
DeferredLogger& deferred_logger) const
{ {
auto pressure_loss = getVfpBhpAdjustment(bhp + dp, thp_limit); auto pressure_loss = getVfpBhpAdjustment(bhp + dp, thp_limit);
auto thp = thp_func(bhp, pressure_loss); auto thp = thp_func(bhp, pressure_loss);
const double tolerance = 1e-5 * unit::barsa; const Scalar tolerance = 1e-5 * unit::barsa;
bool do_iterate = true; bool do_iterate = true;
int it = 1; int it = 1;
int max_iterations = 50; int max_iterations = 50;
while(do_iterate) { while (do_iterate) {
if (it > max_iterations) { if (it > max_iterations) {
break; break;
} }
double thp_prev = thp; Scalar thp_prev = thp;
pressure_loss = getVfpBhpAdjustment(bhp + dp - pressure_loss, thp_prev); pressure_loss = getVfpBhpAdjustment(bhp + dp - pressure_loss, thp_prev);
thp = thp_func(bhp, pressure_loss); thp = thp_func(bhp, pressure_loss);
auto error = std::fabs(thp-thp_prev); auto error = std::fabs(thp - thp_prev);
if (extraThpFromBhpOutput) { if (extraThpFromBhpOutput) {
const std::string msg = fmt::format( const std::string msg = fmt::format(
"findThpFromBhpIteratively(): iteration {}, thp = {}, bhp = {}, " "findThpFromBhpIteratively(): iteration {}, thp = {}, bhp = {}, "
"pressure_loss = {}, error = {}", it, thp, bhp+dp-pressure_loss, pressure_loss, error); "pressure_loss = {}, error = {}", it, thp, bhp+dp-pressure_loss, pressure_loss, error);
deferred_logger.debug(msg); deferred_logger.debug(msg);
} }
if (std::fabs(thp-thp_prev) < tolerance) { if (std::fabs(thp - thp_prev) < tolerance) {
break; break;
} }
it++; it++;
@ -188,14 +193,15 @@ findThpFromBhpIteratively(
return thp; return thp;
} }
std::optional<double> template<class Scalar>
WellBhpThpCalculator:: std::optional<Scalar>
computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates, WellBhpThpCalculator<Scalar>::
computeBhpAtThpLimitProd(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state, const SummaryState& summary_state,
const double maxPerfPress, const Scalar maxPerfPress,
const double rho, const Scalar rho,
const double alq_value, const Scalar alq_value,
const double thp_limit, const Scalar thp_limit,
DeferredLogger& deferred_logger) const DeferredLogger& deferred_logger) const
{ {
// Given a VFP function returning bhp as a function of phase // Given a VFP function returning bhp as a function of phase
@ -221,34 +227,43 @@ computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>&
// Make the fbhp() function. // Make the fbhp() function.
const auto& controls = well_.wellEcl().productionControls(summary_state); const auto& controls = well_.wellEcl().productionControls(summary_state);
const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number); const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const double vfp_ref_depth = table.getDatumDepth(); const Scalar vfp_ref_depth = table.getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); const Scalar 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) { auto fbhp = [this, &controls, thp_limit, dp, alq_value](const std::vector<Scalar>& rates) {
assert(rates.size() == 3); assert(rates.size() == 3);
const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell()); const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number,
const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell()); well_.indexOfWell());
const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number,
well_.indexOfWell());
const bool use_vfpexp = well_.useVfpExplicit(); const bool use_vfpexp = well_.useVfpExplicit();
const double bhp = well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number, const Scalar bhp = well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number,
rates[Water], rates[Water],
rates[Oil], rates[Oil],
rates[Gas], rates[Gas],
thp_limit, thp_limit,
alq_value, alq_value,
wfr, wfr,
gfr, gfr,
use_vfpexp); use_vfpexp);
return bhp - dp + getVfpBhpAdjustment(bhp, thp_limit); return bhp - dp + getVfpBhpAdjustment(bhp, thp_limit);
}; };
// Make the flo() function. // Make the flo() function.
auto flo = [&table](const std::vector<double>& rates) { auto flo = [&table](const std::vector<Scalar>& rates) {
return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]); return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]);
}; };
// Find the bhp-point where production becomes nonzero. // Find the bhp-point where production becomes nonzero.
auto fflo = [&flo, &frates](double bhp) { return flo(frates(bhp)); }; auto fflo = [&flo, &frates](Scalar bhp) { return flo(frates(bhp)); };
auto bhp_max = this->bhpMax(fflo, controls.bhp_limit, maxPerfPress, table.getFloAxis().front(), deferred_logger); 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 // could not solve for the bhp-point, we could not continue to find the bhp
if (!bhp_max.has_value()) { if (!bhp_max.has_value()) {
@ -257,16 +272,17 @@ computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>&
"find bhp-point where production becomes non-zero for well " + well_.name()); "find bhp-point where production becomes non-zero for well " + well_.name());
return std::nullopt; return std::nullopt;
} }
const std::array<double, 2> range {controls.bhp_limit, *bhp_max}; const std::array<Scalar, 2> range {controls.bhp_limit, *bhp_max};
return this->computeBhpAtThpLimit(frates, fbhp, range, deferred_logger); return this->computeBhpAtThpLimit(frates, fbhp, range, deferred_logger);
} }
std::optional<double> template<class Scalar>
WellBhpThpCalculator:: std::optional<Scalar>
computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates, WellBhpThpCalculator<Scalar>::
computeBhpAtThpLimitInj(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state, const SummaryState& summary_state,
const double rho, const Scalar rho,
const double flo_rel_tol, const Scalar flo_rel_tol,
const int max_iteration, const int max_iteration,
const bool throwOnError, const bool throwOnError,
DeferredLogger& deferred_logger) const DeferredLogger& deferred_logger) const
@ -282,13 +298,15 @@ computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>&
} }
} }
void WellBhpThpCalculator::updateThp(const double rho, template<class Scalar>
const bool stop_or_zero_rate_target, void WellBhpThpCalculator<Scalar>::
const std::function<double()>& alq_value, updateThp(const Scalar rho,
const std::array<unsigned,3>& active, const bool stop_or_zero_rate_target,
WellState<double>& well_state, const std::function<Scalar()>& alq_value,
const SummaryState& summary_state, const std::array<unsigned,3>& active,
DeferredLogger& deferred_logger) const WellState<Scalar>& well_state,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const
{ {
static constexpr int Gas = BlackoilPhases::Vapour; static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Oil = BlackoilPhases::Liquid; static constexpr int Oil = BlackoilPhases::Liquid;
@ -308,7 +326,7 @@ void WellBhpThpCalculator::updateThp(const double rho,
} }
// the well is under other control types, we calculate the thp based on bhp and rates // the well is under other control types, we calculate the thp based on bhp and rates
std::vector<double> rates(3, 0.0); std::vector<Scalar> rates(3, 0.0);
const PhaseUsage& pu = well_.phaseUsage(); const PhaseUsage& pu = well_.phaseUsage();
if (active[Water]) { if (active[Water]) {
@ -320,18 +338,19 @@ void WellBhpThpCalculator::updateThp(const double rho,
if (active[Gas]) { if (active[Gas]) {
rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ]; rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ];
} }
const std::optional<double> alq = this->well_.isProducer() ? std::optional<double>(alq_value()) : std::nullopt; const std::optional<Scalar> alq = this->well_.isProducer() ? std::optional<Scalar>(alq_value()) : std::nullopt;
const double thp_limit = well_.getTHPConstraint(summary_state); const Scalar thp_limit = well_.getTHPConstraint(summary_state);
ws.thp = this->calculateThpFromBhp(rates, ws.bhp, rho, alq, thp_limit, deferred_logger); ws.thp = this->calculateThpFromBhp(rates, ws.bhp, rho, alq, thp_limit, deferred_logger);
} }
template<class Scalar>
template<class EvalWell> template<class EvalWell>
EvalWell WellBhpThpCalculator:: EvalWell WellBhpThpCalculator<Scalar>::
calculateBhpFromThp(const WellState<double>& well_state, calculateBhpFromThp(const WellState<Scalar>& well_state,
const std::vector<EvalWell>& rates, const std::vector<EvalWell>& rates,
const Well& well, const Well& well,
const SummaryState& summaryState, const SummaryState& summaryState,
const double rho, const Scalar rho,
DeferredLogger& deferred_logger) const DeferredLogger& deferred_logger) const
{ {
// TODO: when well is under THP control, the BHP is dependent on the rates, // TODO: when well is under THP control, the BHP is dependent on the rates,
@ -349,8 +368,8 @@ calculateBhpFromThp(const WellState<double>& well_state,
const EvalWell aqua = rates[Water]; const EvalWell aqua = rates[Water];
const EvalWell liquid = rates[Oil]; const EvalWell liquid = rates[Oil];
const EvalWell vapour = rates[Gas]; const EvalWell vapour = rates[Gas];
const double thp_limit = well_.getTHPConstraint(summaryState); const Scalar thp_limit = well_.getTHPConstraint(summaryState);
double vfp_ref_depth; Scalar vfp_ref_depth;
EvalWell bhp_tab; EvalWell bhp_tab;
if (well_.isInjector() ) if (well_.isInjector() )
{ {
@ -375,57 +394,62 @@ calculateBhpFromThp(const WellState<double>& well_state,
else { else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + well_.name(), deferred_logger); OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + well_.name(), deferred_logger);
} }
double bhp_tab_double_value; Scalar bhp_tab_double_value;
if constexpr (std::is_same_v<EvalWell, double>) { if constexpr (std::is_same_v<EvalWell, Scalar>) {
bhp_tab_double_value = bhp_tab; bhp_tab_double_value = bhp_tab;
} }
else { // EvalWell and bhp_tab is of type DenseAd::Evaluation<double,...,...> else { // EvalWell and bhp_tab is of type DenseAd::Evaluation<double,...,...>
bhp_tab_double_value = bhp_tab.value(); bhp_tab_double_value = bhp_tab.value();
} }
const auto bhp_adjustment = getVfpBhpAdjustment(bhp_tab_double_value, thp_limit); const auto bhp_adjustment = getVfpBhpAdjustment(bhp_tab_double_value, thp_limit);
const double dp_hydro = wellhelpers::computeHydrostaticCorrection( const Scalar dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(),
well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); vfp_ref_depth,
rho,
well_.gravity());
return bhp_tab - dp_hydro + bhp_adjustment; return bhp_tab - dp_hydro + bhp_adjustment;
} }
double template<class Scalar>
WellBhpThpCalculator:: Scalar WellBhpThpCalculator<Scalar>::
calculateMinimumBhpFromThp(const WellState<double>& well_state, calculateMinimumBhpFromThp(const WellState<Scalar>& well_state,
const Well& well, const Well& well,
const SummaryState& summaryState, const SummaryState& summaryState,
const double rho) const const Scalar rho) const
{ {
assert(well_.isProducer()); // only producers can go here for now assert(well_.isProducer()); // only producers can go here for now
const double thp_limit = well_.getTHPConstraint(summaryState); const Scalar thp_limit = well_.getTHPConstraint(summaryState);
const auto& controls = well.productionControls(summaryState); const auto& controls = well.productionControls(summaryState);
const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell()); 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 auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell());
const double bhp_min = well_.vfpProperties()->getProd()->minimumBHP(controls.vfp_table_number, const Scalar bhp_min = well_.vfpProperties()->getProd()->minimumBHP(controls.vfp_table_number,
thp_limit, wfr, gfr, thp_limit, wfr, gfr,
well_.getALQ(well_state)); well_.getALQ(well_state));
const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth(); const Scalar vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const auto bhp_adjustment = getVfpBhpAdjustment(bhp_min, thp_limit); const auto bhp_adjustment = getVfpBhpAdjustment(bhp_min, thp_limit);
const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, const Scalar dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth,
rho, well_.gravity()); rho, well_.gravity());
return bhp_min - dp_hydro + bhp_adjustment; return bhp_min - dp_hydro + bhp_adjustment;
} }
double WellBhpThpCalculator::getVfpBhpAdjustment(const double bhp_tab, const double thp_limit) const template<class Scalar>
Scalar WellBhpThpCalculator<Scalar>::
getVfpBhpAdjustment(const Scalar bhp_tab, const Scalar thp_limit) const
{ {
return well_.wellEcl().getWVFPDP().getPressureLoss(bhp_tab, thp_limit); return well_.wellEcl().getWVFPDP().getPressureLoss(bhp_tab, thp_limit);
} }
template<class Scalar>
template<class ErrorPolicy> template<class ErrorPolicy>
std::optional<double> std::optional<Scalar>
WellBhpThpCalculator:: WellBhpThpCalculator<Scalar>::
computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double)>& frates, computeBhpAtThpLimitInjImpl(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state, const SummaryState& summary_state,
const double rho, const Scalar rho,
const double flo_rel_tol, const Scalar flo_rel_tol,
const int max_iteration, const int max_iteration,
DeferredLogger& deferred_logger) const DeferredLogger& deferred_logger) const
{ {
@ -475,12 +499,12 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
// Make the fbhp() function. // Make the fbhp() function.
const auto& controls = well_.wellEcl().injectionControls(summary_state); const auto& controls = well_.wellEcl().injectionControls(summary_state);
const auto& table = well_.vfpProperties()->getInj()->getTable(controls.vfp_table_number); const auto& table = well_.vfpProperties()->getInj()->getTable(controls.vfp_table_number);
const double vfp_ref_depth = table.getDatumDepth(); const Scalar vfp_ref_depth = table.getDatumDepth();
const double thp_limit = this->getTHPConstraint(summary_state); const Scalar thp_limit = this->getTHPConstraint(summary_state);
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(),
vfp_ref_depth, vfp_ref_depth,
rho, well_.gravity()); rho, well_.gravity());
auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& rates) { auto fbhp = [this, &controls, thp_limit, dp](const std::vector<Scalar>& rates) {
assert(rates.size() == 3); assert(rates.size() == 3);
const auto bhp = well_.vfpProperties()->getInj() const auto bhp = well_.vfpProperties()->getInj()
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit); ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit);
@ -488,25 +512,25 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
}; };
// Make the flo() function. // Make the flo() function.
auto flo = [&table](const std::vector<double>& rates) { auto flo = [&table](const std::vector<Scalar>& rates) {
return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]); return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]);
}; };
// Get the flo samples, add extra samples at low rates and bhp // Get the flo samples, add extra samples at low rates and bhp
// limit point if necessary. // limit point if necessary.
std::vector<double> flo_samples = table.getFloAxis(); std::vector<Scalar> flo_samples = table.getFloAxis();
if (flo_samples[0] > 0.0) { if (flo_samples[0] > 0.0) {
const double f0 = flo_samples[0]; const Scalar f0 = flo_samples[0];
flo_samples.insert(flo_samples.begin(), { f0/20.0, f0/10.0, f0/5.0, f0/2.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)); const Scalar flo_bhp_limit = flo(frates(controls.bhp_limit));
if (flo_samples.back() < flo_bhp_limit) { if (flo_samples.back() < flo_bhp_limit) {
flo_samples.push_back(flo_bhp_limit); flo_samples.push_back(flo_bhp_limit);
} }
// Find bhp values for inflow relation corresponding to flo samples. // Find bhp values for inflow relation corresponding to flo samples.
std::vector<double> bhp_samples; std::vector<Scalar> bhp_samples;
for (double flo_sample : flo_samples) { for (Scalar flo_sample : flo_samples) {
if (flo_sample > flo_bhp_limit) { if (flo_sample > flo_bhp_limit) {
// We would have to go over the bhp limit to obtain a // We would have to go over the bhp limit to obtain a
// flow of this magnitude. We associate all such flows // flow of this magnitude. We associate all such flows
@ -516,16 +540,16 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
bhp_samples.push_back(controls.bhp_limit); bhp_samples.push_back(controls.bhp_limit);
break; break;
} }
auto eq = [&flo, &frates, flo_sample](double bhp) { auto eq = [&flo, &frates, flo_sample](Scalar bhp) {
return flo(frates(bhp)) - flo_sample; return flo(frates(bhp)) - flo_sample;
}; };
// TODO: replace hardcoded low/high limits. // TODO: replace hardcoded low/high limits.
const double low = 10.0 * unit::barsa; const Scalar low = 10.0 * unit::barsa;
const double high = 800.0 * unit::barsa; const Scalar high = 800.0 * unit::barsa;
const double flo_tolerance = flo_rel_tol * std::fabs(flo_samples.back()); const Scalar flo_tolerance = flo_rel_tol * std::fabs(flo_samples.back());
int iteration = 0; int iteration = 0;
try { try {
const double solved_bhp = RegulaFalsiBisection<ErrorPolicy>:: const Scalar solved_bhp = RegulaFalsiBisection<ErrorPolicy>::
solve(eq, low, high, max_iteration, flo_tolerance, iteration); solve(eq, low, high, max_iteration, flo_tolerance, iteration);
bhp_samples.push_back(solved_bhp); bhp_samples.push_back(solved_bhp);
} }
@ -539,7 +563,7 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
// Find bhp values for VFP relation corresponding to flo samples. // 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() const int num_samples = bhp_samples.size(); // Note that this can be smaller than flo_samples.size()
std::vector<double> fbhp_samples(num_samples); std::vector<Scalar> fbhp_samples(num_samples);
for (int ii = 0; ii < num_samples; ++ii) { for (int ii = 0; ii < num_samples; ++ii) {
fbhp_samples[ii] = fbhp(frates(bhp_samples[ii])); fbhp_samples[ii] = fbhp(frates(bhp_samples[ii]));
} }
@ -564,8 +588,8 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
// We only look at the valid // We only look at the valid
int sign_change_index = -1; int sign_change_index = -1;
for (int ii = 0; ii < num_samples - 1; ++ii) { for (int ii = 0; ii < num_samples - 1; ++ii) {
const double curr = fbhp_samples[ii] - bhp_samples[ii]; const Scalar curr = fbhp_samples[ii] - bhp_samples[ii];
const double next = fbhp_samples[ii + 1] - bhp_samples[ii + 1]; const Scalar next = fbhp_samples[ii + 1] - bhp_samples[ii + 1];
if (curr * next < 0.0) { if (curr * next < 0.0) {
// Sign change in the [ii, ii + 1] interval. // Sign change in the [ii, ii + 1] interval.
sign_change_index = ii; // May overwrite, thereby choosing the highest-flo solution. sign_change_index = ii; // May overwrite, thereby choosing the highest-flo solution.
@ -578,13 +602,13 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
} }
// Solve for the proper solution in the given interval. // Solve for the proper solution in the given interval.
auto eq = [&fbhp, &frates](double bhp) { auto eq = [&fbhp, &frates](Scalar bhp) {
return fbhp(frates(bhp)) - bhp; return fbhp(frates(bhp)) - bhp;
}; };
// TODO: replace hardcoded low/high limits. // TODO: replace hardcoded low/high limits.
const double low = bhp_samples[sign_change_index + 1]; const Scalar low = bhp_samples[sign_change_index + 1];
const double high = bhp_samples[sign_change_index]; const Scalar high = bhp_samples[sign_change_index];
const double bhp_tolerance = 0.01 * unit::barsa; const Scalar bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0; int iteration = 0;
if (low == high) { if (low == high) {
// We are in the high flow regime where the bhp_samples // We are in the high flow regime where the bhp_samples
@ -595,7 +619,7 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
return std::nullopt; return std::nullopt;
} }
try { try {
const double solved_bhp = RegulaFalsiBisection<ErrorPolicy>:: const Scalar solved_bhp = RegulaFalsiBisection<ErrorPolicy>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration); solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
if constexpr (extraBhpAtThpLimitOutput) { if constexpr (extraBhpAtThpLimitOutput) {
OpmLog::debug("***** " + well_.name() + " solved_bhp = " + std::to_string(solved_bhp) OpmLog::debug("***** " + well_.name() + " solved_bhp = " + std::to_string(solved_bhp)
@ -610,20 +634,21 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
} }
} }
std::optional<double> template<class Scalar>
WellBhpThpCalculator:: std::optional<Scalar>
bhpMax(const std::function<double(const double)>& fflo, WellBhpThpCalculator<Scalar>::
const double bhp_limit, bhpMax(const std::function<Scalar(const Scalar)>& fflo,
const double maxPerfPress, const Scalar bhp_limit,
const double vfp_flo_front, const Scalar maxPerfPress,
const Scalar vfp_flo_front,
DeferredLogger& deferred_logger) const DeferredLogger& deferred_logger) const
{ {
// Find the bhp-point where production becomes nonzero. // Find the bhp-point where production becomes nonzero.
double bhp_max = 0.0; Scalar bhp_max = 0.0;
double low = bhp_limit; Scalar low = bhp_limit;
double high = maxPerfPress + 1.0 * unit::barsa; Scalar high = maxPerfPress + 1.0 * unit::barsa;
double f_low = fflo(low); Scalar f_low = fflo(low);
double f_high = fflo(high); Scalar f_high = fflo(high);
if constexpr (extraBhpAtThpLimitOutput) { if constexpr (extraBhpAtThpLimitOutput) {
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + well_.name() + deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + well_.name() +
" low = " + std::to_string(low) + " low = " + std::to_string(low) +
@ -633,7 +658,7 @@ bhpMax(const std::function<double(const double)>& fflo,
} }
int adjustments = 0; int adjustments = 0;
const int max_adjustments = 10; const int max_adjustments = 10;
const double adjust_amount = 5.0 * unit::barsa; const Scalar adjust_amount = 5.0 * unit::barsa;
while (f_low * f_high > 0.0 && adjustments < max_adjustments) { while (f_low * f_high > 0.0 && adjustments < max_adjustments) {
// Same sign, adjust high to see if we can flip it. // Same sign, adjust high to see if we can flip it.
high += adjust_amount; high += adjust_amount;
@ -656,12 +681,12 @@ bhpMax(const std::function<double(const double)>& fflo,
} else { } else {
// Bisect to find a bhp point where we produce, but // Bisect to find a bhp point where we produce, but
// not a large amount ('eps' below). // not a large amount ('eps' below).
const double eps = 0.1 * std::fabs(vfp_flo_front); const Scalar eps = 0.1 * std::fabs(vfp_flo_front);
const int maxit = 50; const int maxit = 50;
int it = 0; int it = 0;
while (std::fabs(f_low) > eps && it < maxit) { while (std::fabs(f_low) > eps && it < maxit) {
const double curr = 0.5*(low + high); const Scalar curr = 0.5*(low + high);
const double f_curr = fflo(curr); const Scalar f_curr = fflo(curr);
if (f_curr * f_low > 0.0) { if (f_curr * f_low > 0.0) {
low = curr; low = curr;
f_low = f_curr; f_low = f_curr;
@ -690,11 +715,12 @@ bhpMax(const std::function<double(const double)>& fflo,
return bhp_max; return bhp_max;
} }
std::optional<double> template<class Scalar>
WellBhpThpCalculator:: std::optional<Scalar>
computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& frates, WellBhpThpCalculator<Scalar>::
const std::function<double(const std::vector<double>)>& fbhp, computeBhpAtThpLimit(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const std::array<double, 2>& range, const std::function<Scalar(const std::vector<Scalar>)>& fbhp,
const std::array<Scalar, 2>& range,
DeferredLogger& deferred_logger) const DeferredLogger& deferred_logger) const
{ {
// Given a VFP function returning bhp as a function of phase // Given a VFP function returning bhp as a function of phase
@ -714,15 +740,16 @@ computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& fra
// highest rate) should be returned. // highest rate) should be returned.
// Define the equation we want to solve. // Define the equation we want to solve.
auto eq = [&fbhp, &frates](double bhp) { auto eq = [&fbhp, &frates](Scalar bhp) {
return fbhp(frates(bhp)) - bhp; return fbhp(frates(bhp)) - bhp;
}; };
// Find appropriate brackets for the solution. // Find appropriate brackets for the solution.
std::optional<double> approximate_solution; std::optional<Scalar> approximate_solution;
double low, high; Scalar low, high;
// trying to use bisect way to locate a bracket // trying to use bisect way to locate a bracket
bool finding_bracket = this->bisectBracket(eq, range, low, high, approximate_solution, deferred_logger); 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 // based on the origional design, if an approximate solution is suggested, we use this value directly
// in the long run, we might change it // in the long run, we might change it
@ -744,10 +771,10 @@ computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& fra
// Solve for the proper solution in the given interval. // Solve for the proper solution in the given interval.
const int max_iteration = 100; const int max_iteration = 100;
const double bhp_tolerance = 0.01 * unit::barsa; const Scalar bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0; int iteration = 0;
try { try {
const double solved_bhp = RegulaFalsiBisection<ThrowOnError>:: const Scalar solved_bhp = RegulaFalsiBisection<ThrowOnError>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration); solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
return solved_bhp; return solved_bhp;
} }
@ -758,21 +785,21 @@ computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& fra
} }
} }
bool template<class Scalar>
WellBhpThpCalculator:: bool WellBhpThpCalculator<Scalar>::
bisectBracket(const std::function<double(const double)>& eq, bisectBracket(const std::function<Scalar(const Scalar)>& eq,
const std::array<double, 2>& range, const std::array<Scalar, 2>& range,
double& low, double& high, Scalar& low, Scalar& high,
std::optional<double>& approximate_solution, std::optional<Scalar>& approximate_solution,
DeferredLogger& deferred_logger) const DeferredLogger& deferred_logger) const
{ {
bool finding_bracket = false; bool finding_bracket = false;
low = range[0]; low = range[0];
high = range[1]; high = range[1];
double eq_high = eq(high); Scalar eq_high = eq(high);
double eq_low = eq(low); Scalar eq_low = eq(low);
const double eq_bhplimit = eq_low; const Scalar eq_bhplimit = eq_low;
if constexpr (extraBhpAtThpLimitOutput) { if constexpr (extraBhpAtThpLimitOutput) {
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + well_.name() + deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + well_.name() +
" low = " + std::to_string(low) + " low = " + std::to_string(low) +
@ -783,12 +810,12 @@ bisectBracket(const std::function<double(const double)>& eq,
if (eq_low * eq_high > 0.0) { if (eq_low * eq_high > 0.0) {
// Failed to bracket the zero. // Failed to bracket the zero.
// If this is due to having two solutions, bisect until bracketed. // If this is due to having two solutions, bisect until bracketed.
double abs_low = std::fabs(eq_low); Scalar abs_low = std::fabs(eq_low);
double abs_high = std::fabs(eq_high); Scalar abs_high = std::fabs(eq_high);
int bracket_attempts = 0; int bracket_attempts = 0;
const int max_bracket_attempts = 20; const int max_bracket_attempts = 20;
double interval = high - low; Scalar interval = high - low;
const double min_interval = 1.0 * unit::barsa; const Scalar min_interval = 1.0 * unit::barsa;
while (eq_low * eq_high > 0.0 && bracket_attempts < max_bracket_attempts && interval > min_interval) { while (eq_low * eq_high > 0.0 && bracket_attempts < max_bracket_attempts && interval > min_interval) {
if (abs_high < abs_low) { if (abs_high < abs_low) {
low = 0.5 * (low + high); low = 0.5 * (low + high);
@ -815,7 +842,7 @@ bisectBracket(const std::function<double(const double)>& eq,
} }
} else { // eq_low * eq_high > 0.0 } else { // eq_low * eq_high > 0.0
// Still failed bracketing! // Still failed bracketing!
const double limit = 0.1 * unit::barsa; const Scalar limit = 0.1 * unit::barsa;
if (std::min(abs_low, abs_high) < limit) { if (std::min(abs_low, abs_high) < limit) {
// Return the least bad solution if less off than 0.1 bar. // Return the least bad solution if less off than 0.1 bar.
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE", deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE",
@ -834,20 +861,20 @@ bisectBracket(const std::function<double(const double)>& eq,
return finding_bracket; return finding_bracket;
} }
bool template<class Scalar>
WellBhpThpCalculator:: bool WellBhpThpCalculator<Scalar>::
bruteForceBracket(const std::function<double(const double)>& eq, bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
const std::array<double, 2>& range, const std::array<Scalar, 2>& range,
double& low, double& high, Scalar& low, Scalar& high,
DeferredLogger& deferred_logger) DeferredLogger& deferred_logger)
{ {
bool bracket_found = false; bool bracket_found = false;
low = range[0]; low = range[0];
high = range[1]; high = range[1];
const int sample_number = 200; const int sample_number = 200;
const double interval = (high - low) / sample_number; const Scalar interval = (high - low) / sample_number;
double eq_low = eq(low); Scalar eq_low = eq(low);
double eq_high = 0.0; Scalar eq_high = 0.0;
for (int i = 0; i < sample_number + 1; ++i) { for (int i = 0; i < sample_number + 1; ++i) {
high = range[0] + interval * i; high = range[0] + interval * i;
eq_high = eq(high); eq_high = eq(high);
@ -866,10 +893,11 @@ bruteForceBracket(const std::function<double(const double)>& eq,
return bracket_found; return bracket_found;
} }
bool WellBhpThpCalculator:: template<class Scalar>
isStableSolution(const WellState<double>& well_state, bool WellBhpThpCalculator<Scalar>::
isStableSolution(const WellState<Scalar>& well_state,
const Well& well, const Well& well,
const std::vector<double>& rates, const std::vector<Scalar>& rates,
const SummaryState& summaryState) const const SummaryState& summaryState) const
{ {
assert(int(rates.size()) == 3); // the vfp related only supports three phases now. assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
@ -879,10 +907,10 @@ isStableSolution(const WellState<double>& well_state,
static constexpr int Oil = BlackoilPhases::Liquid; static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua; static constexpr int Water = BlackoilPhases::Aqua;
const double aqua = rates[Water]; const Scalar aqua = rates[Water];
const double liquid = rates[Oil]; const Scalar liquid = rates[Oil];
const double vapour = rates[Gas]; const Scalar vapour = rates[Gas];
const double thp = well_.getTHPConstraint(summaryState); const Scalar thp = well_.getTHPConstraint(summaryState);
const auto& controls = well.productionControls(summaryState); const auto& controls = well.productionControls(summaryState);
const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell()); const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
@ -897,27 +925,28 @@ isStableSolution(const WellState<double>& well_state,
return true; return true;
} else { // maybe check if ipr is available } else { // maybe check if ipr is available
const auto ipr = getFloIPR(well_state, well, summaryState); const auto ipr = getFloIPR(well_state, well, summaryState);
return bhp.dflo + 1/ipr.second >= 0; return bhp.dflo + 1.0 / ipr.second >= 0;
} }
} }
std::optional<double> WellBhpThpCalculator:: template<class Scalar>
estimateStableBhp(const WellState<double>& well_state, std::optional<Scalar> WellBhpThpCalculator<Scalar>::
estimateStableBhp(const WellState<Scalar>& well_state,
const Well& well, const Well& well,
const std::vector<double>& rates, const std::vector<Scalar>& rates,
const double rho, const Scalar rho,
const SummaryState& summaryState) const const SummaryState& summaryState) const
{ {
// Given a *converged* well_state with ipr, estimate bhp of the stable solution // Given a *converged* well_state with ipr, estimate bhp of the stable solution
const auto& controls = well.productionControls(summaryState); const auto& controls = well.productionControls(summaryState);
const double thp = well_.getTHPConstraint(summaryState); const Scalar thp = well_.getTHPConstraint(summaryState);
const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number); const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const double aqua = rates[BlackoilPhases::Aqua]; const Scalar aqua = rates[BlackoilPhases::Aqua];
const double liquid = rates[BlackoilPhases::Liquid]; const Scalar liquid = rates[BlackoilPhases::Liquid];
const double vapour = rates[BlackoilPhases::Vapour]; const Scalar vapour = rates[BlackoilPhases::Vapour];
double flo = detail::getFlo(table, aqua, liquid, vapour); Scalar flo = detail::getFlo(table, aqua, liquid, vapour);
double wfr, gfr; Scalar wfr, gfr;
if (well_.useVfpExplicit() || -flo < table.getFloAxis().front()) { if (well_.useVfpExplicit() || -flo < table.getFloAxis().front()) {
wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell()); wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell()); gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell());
@ -928,11 +957,11 @@ estimateStableBhp(const WellState<double>& well_state,
auto ipr = getFloIPR(well_state, well, summaryState); auto ipr = getFloIPR(well_state, well, summaryState);
const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth(); const Scalar vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, const Scalar dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth,
rho, well_.gravity()); rho, well_.gravity());
auto bhp_adjusted = [this, &thp, &dp_hydro](const double bhp) { auto bhp_adjusted = [this, &thp, &dp_hydro](const Scalar bhp) {
return bhp - dp_hydro + getVfpBhpAdjustment(bhp, thp); return bhp - dp_hydro + getVfpBhpAdjustment(bhp, thp);
}; };
const auto retval = detail::intersectWithIPR(table, thp, wfr, gfr, well_.getALQ(well_state), ipr.first, ipr.second, bhp_adjusted); const auto retval = detail::intersectWithIPR(table, thp, wfr, gfr, well_.getALQ(well_state), ipr.first, ipr.second, bhp_adjusted);
@ -944,8 +973,9 @@ estimateStableBhp(const WellState<double>& well_state,
} }
} }
std::pair<double, double> WellBhpThpCalculator:: template<class Scalar>
getFloIPR(const WellState<double>& well_state, std::pair<Scalar, Scalar> WellBhpThpCalculator<Scalar>::
getFloIPR(const WellState<Scalar>& well_state,
const Well& well, const Well& well,
const SummaryState& summary_state) const const SummaryState& summary_state) const
{ {
@ -954,21 +984,23 @@ getFloIPR(const WellState<double>& well_state,
const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number); const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const auto& pu = well_.phaseUsage(); const auto& pu = well_.phaseUsage();
const auto& ipr_a = well_state.well(well_.indexOfWell()).implicit_ipr_a; const auto& ipr_a = well_state.well(well_.indexOfWell()).implicit_ipr_a;
const double& aqua_a = pu.phase_used[BlackoilPhases::Aqua]? ipr_a[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0; const Scalar& aqua_a = pu.phase_used[BlackoilPhases::Aqua]? ipr_a[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0;
const double& liquid_a = pu.phase_used[BlackoilPhases::Liquid]? ipr_a[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0; const Scalar& liquid_a = pu.phase_used[BlackoilPhases::Liquid]? ipr_a[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0;
const double& vapour_a = pu.phase_used[BlackoilPhases::Vapour]? ipr_a[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0; const Scalar& vapour_a = pu.phase_used[BlackoilPhases::Vapour]? ipr_a[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0;
const auto& ipr_b = well_state.well(well_.indexOfWell()).implicit_ipr_b; const auto& ipr_b = well_state.well(well_.indexOfWell()).implicit_ipr_b;
const double& aqua_b = pu.phase_used[BlackoilPhases::Aqua]? ipr_b[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0; const Scalar& aqua_b = pu.phase_used[BlackoilPhases::Aqua]? ipr_b[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0;
const double& liquid_b = pu.phase_used[BlackoilPhases::Liquid]? ipr_b[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0; const Scalar& liquid_b = pu.phase_used[BlackoilPhases::Liquid]? ipr_b[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0;
const double& vapour_b = pu.phase_used[BlackoilPhases::Vapour]? ipr_b[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0; const Scalar& vapour_b = pu.phase_used[BlackoilPhases::Vapour]? ipr_b[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0;
// The getFlo helper is indended to pick one or add two of the phase rates (depending on FLO-type), // The getFlo helper is indended to pick one or add two of the phase rates (depending on FLO-type),
// but we can equally use it to pick/add the corresponding ipr_a, ipr_b // but we can equally use it to pick/add the corresponding ipr_a, ipr_b
return std::make_pair(detail::getFlo(table, aqua_a, liquid_a, vapour_a), return std::make_pair(detail::getFlo(table, aqua_a, liquid_a, vapour_a),
detail::getFlo(table, aqua_b, liquid_b, vapour_b)); detail::getFlo(table, aqua_b, liquid_b, vapour_b));
} }
template class WellBhpThpCalculator<double>;
#define INSTANCE(...) \ #define INSTANCE(...) \
template __VA_ARGS__ WellBhpThpCalculator:: \ template __VA_ARGS__ WellBhpThpCalculator<double>:: \
calculateBhpFromThp<__VA_ARGS__>(const WellState<double>&, \ calculateBhpFromThp<__VA_ARGS__>(const WellState<double>&, \
const std::vector<__VA_ARGS__>&, \ const std::vector<__VA_ARGS__>&, \
const Well&, \ const Well&, \
@ -993,4 +1025,5 @@ INSTANCE(DenseAd::Evaluation<double,-1,8u>)
INSTANCE(DenseAd::Evaluation<double,-1,9u>) INSTANCE(DenseAd::Evaluation<double,-1,9u>)
INSTANCE(DenseAd::Evaluation<double,-1,10u>) INSTANCE(DenseAd::Evaluation<double,-1,10u>)
INSTANCE(DenseAd::Evaluation<double,-1,11u>) INSTANCE(DenseAd::Evaluation<double,-1,11u>)
} // namespace Opm } // namespace Opm

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

@ -26,11 +26,9 @@
#include <functional> #include <functional>
#include <optional> #include <optional>
#include <string>
#include <vector> #include <vector>
namespace Opm namespace Opm {
{
class DeferredLogger; class DeferredLogger;
class SummaryState; class SummaryState;
@ -39,136 +37,137 @@ template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState; template<class Scalar> class WellState;
//! \brief Class for computing BHP limits. //! \brief Class for computing BHP limits.
template<class Scalar>
class WellBhpThpCalculator { class WellBhpThpCalculator {
public: public:
//! \brief Constructor sets reference to well. //! \brief Constructor sets reference to well.
WellBhpThpCalculator(const WellInterfaceGeneric<double>& well) : well_(well) {} WellBhpThpCalculator(const WellInterfaceGeneric<Scalar>& well) : well_(well) {}
//! \brief Checks if well has THP constraints. //! \brief Checks if well has THP constraints.
bool wellHasTHPConstraints(const SummaryState& summaryState) const; bool wellHasTHPConstraints(const SummaryState& summaryState) const;
//! \brief Get THP constraint for well. //! \brief Get THP constraint for well.
double getTHPConstraint(const SummaryState& summaryState) const; Scalar getTHPConstraint(const SummaryState& summaryState) const;
//! \brief Obtain the most strict BHP from BHP limits. //! \brief Obtain the most strict BHP from BHP limits.
double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const; Scalar mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
//! \brief Calculates THP from BHP. //! \brief Calculates THP from BHP.
double calculateThpFromBhp(const std::vector<double>& rates, Scalar calculateThpFromBhp(const std::vector<Scalar>& rates,
const double bhp, const Scalar bhp,
const double rho, const Scalar rho,
const std::optional<double>& alq, const std::optional<Scalar>& alq,
const double thp_limit, const Scalar thp_limit,
DeferredLogger& deferred_logger) const; DeferredLogger& deferred_logger) const;
//! \brief Compute BHP from THP limit for a producer. //! \brief Compute BHP from THP limit for a producer.
std::optional<double> std::optional<Scalar>
computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates, computeBhpAtThpLimitProd(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state, const SummaryState& summary_state,
const double maxPerfPress, const Scalar maxPerfPress,
const double rho, const Scalar rho,
const double alq_value, const Scalar alq_value,
const double thp_limit, const Scalar thp_limit,
DeferredLogger& deferred_logger) const; DeferredLogger& deferred_logger) const;
//! \brief Compute BHP from THP limit for an injector. //! \brief Compute BHP from THP limit for an injector.
std::optional<double> std::optional<Scalar>
computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates, computeBhpAtThpLimitInj(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state, const SummaryState& summary_state,
const double rho, const Scalar rho,
const double flo_rel_tol, const Scalar flo_rel_tol,
const int max_iteration, const int max_iteration,
const bool throwOnError, const bool throwOnError,
DeferredLogger& deferred_logger) const; DeferredLogger& deferred_logger) const;
//! \brief Update THP. //! \brief Update THP.
void updateThp(const double rho, void updateThp(const Scalar rho,
const bool stop_or_zero_rate_target, const bool stop_or_zero_rate_target,
const std::function<double()>& alq_value, const std::function<Scalar()>& alq_value,
const std::array<unsigned,3>& active, const std::array<unsigned,3>& active,
WellState<double>& well_state, WellState<Scalar>& well_state,
const SummaryState& summary_state, const SummaryState& summary_state,
DeferredLogger& deferred_logger) const; DeferredLogger& deferred_logger) const;
template<class EvalWell> template<class EvalWell>
EvalWell calculateBhpFromThp(const WellState<double>& well_state, EvalWell calculateBhpFromThp(const WellState<Scalar>& well_state,
const std::vector<EvalWell>& rates, const std::vector<EvalWell>& rates,
const Well& well, const Well& well,
const SummaryState& summaryState, const SummaryState& summaryState,
const double rho, const Scalar rho,
DeferredLogger& deferred_logger) const; DeferredLogger& deferred_logger) const;
double calculateMinimumBhpFromThp(const WellState<double>& well_state, Scalar calculateMinimumBhpFromThp(const WellState<Scalar>& well_state,
const Well& well, const Well& well,
const SummaryState& summaryState, const SummaryState& summaryState,
const double rho) const; const Scalar rho) const;
bool isStableSolution(const WellState<double>& well_state, bool isStableSolution(const WellState<Scalar>& well_state,
const Well& well, const Well& well,
const std::vector<double>& rates, const std::vector<Scalar>& rates,
const SummaryState& summaryState) const; const SummaryState& summaryState) const;
std::optional<double> std::optional<Scalar>
estimateStableBhp (const WellState<double>& well_state, estimateStableBhp (const WellState<Scalar>& well_state,
const Well& well, const Well& well,
const std::vector<double>& rates, const std::vector<Scalar>& rates,
const double rho, const Scalar rho,
const SummaryState& summaryState) const; const SummaryState& summaryState) const;
std::pair<double, double> std::pair<Scalar, Scalar>
getFloIPR(const WellState<double>& well_state, getFloIPR(const WellState<Scalar>& well_state,
const Well& well, const Well& well,
const SummaryState& summary_state) const; const SummaryState& summary_state) const;
private: private:
//! \brief Compute BHP from THP limit for an injector - implementation. //! \brief Compute BHP from THP limit for an injector - implementation.
template<class ErrorPolicy> template<class ErrorPolicy>
std::optional<double> std::optional<Scalar>
computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double)>& frates, computeBhpAtThpLimitInjImpl(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state, const SummaryState& summary_state,
const double rho, const Scalar rho,
const double flo_rel_tol, const Scalar flo_rel_tol,
const int max_iteration, const int max_iteration,
DeferredLogger& deferred_logger) const; DeferredLogger& deferred_logger) const;
//! \brief Calculate max BHP. //! \brief Calculate max BHP.
std::optional<double> std::optional<Scalar>
bhpMax(const std::function<double(const double)>& fflo, bhpMax(const std::function<Scalar(const Scalar)>& fflo,
const double bhp_limit, const Scalar bhp_limit,
const double maxPerfPress, const Scalar maxPerfPress,
const double vfp_flo_front, const Scalar vfp_flo_front,
DeferredLogger& deferred_logger) const; DeferredLogger& deferred_logger) const;
//! \brief Common code for finding BHP from THP limit for producers/injectors. //! \brief Common code for finding BHP from THP limit for producers/injectors.
std::optional<double> std::optional<Scalar>
computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& frates, computeBhpAtThpLimit(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const std::function<double(const std::vector<double>)>& fbhp, const std::function<Scalar(const std::vector<Scalar>)>& fbhp,
const std::array<double, 2>& range, const std::array<Scalar, 2>& range,
DeferredLogger& deferred_logger) const; DeferredLogger& deferred_logger) const;
//! \brief Get pressure adjustment to the bhp calculated from VFP table //! \brief Get pressure adjustment to the bhp calculated from VFP table
double getVfpBhpAdjustment(const double bph_tab, const double thp_limit) const; Scalar getVfpBhpAdjustment(const Scalar bph_tab, const Scalar thp_limit) const;
//! \brief Find limits using bisection. //! \brief Find limits using bisection.
bool bisectBracket(const std::function<double(const double)>& eq, bool bisectBracket(const std::function<Scalar(const Scalar)>& eq,
const std::array<double, 2>& range, const std::array<Scalar, 2>& range,
double& low, double& high, Scalar& low, Scalar& high,
std::optional<double>& approximate_solution, std::optional<Scalar>& approximate_solution,
DeferredLogger& deferred_logger) const; DeferredLogger& deferred_logger) const;
//! \brief Find limits using brute-force solver. //! \brief Find limits using brute-force solver.
static bool bruteForceBracket(const std::function<double(const double)>& eq, static bool bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
const std::array<double, 2>& range, const std::array<Scalar, 2>& range,
double& low, double& high, Scalar& low, Scalar& high,
DeferredLogger& deferred_logger); DeferredLogger& deferred_logger);
double findThpFromBhpIteratively(const std::function<double(const double, const double)>& thp_func, Scalar findThpFromBhpIteratively(const std::function<Scalar(const Scalar, const Scalar)>& thp_func,
const double bhp, const Scalar bhp,
const double thp_limit, const Scalar thp_limit,
const double dp, const Scalar dp,
DeferredLogger& deferred_logger) const; DeferredLogger& deferred_logger) const;
const WellInterfaceGeneric<double>& well_; //!< Reference to well interface const WellInterfaceGeneric<Scalar>& well_; //!< Reference to well interface
}; };
} }