From 1001d35418e854ed14f778932c3865989945501c Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 9 Feb 2022 10:06:09 +0100 Subject: [PATCH] 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 --- .../wells/MultisegmentWellGeneric.cpp | 279 +---------------- opm/simulators/wells/StandardWellGeneric.cpp | 5 + opm/simulators/wells/StandardWellGeneric.hpp | 1 + opm/simulators/wells/StandardWell_impl.hpp | 9 + opm/simulators/wells/WellInterfaceGeneric.cpp | 287 +++++++++++++++++- opm/simulators/wells/WellInterfaceGeneric.hpp | 24 ++ 6 files changed, 326 insertions(+), 279 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellGeneric.cpp b/opm/simulators/wells/MultisegmentWellGeneric.cpp index 48fa91054..ae98d5a79 100644 --- a/opm/simulators/wells/MultisegmentWellGeneric.cpp +++ b/opm/simulators/wells/MultisegmentWellGeneric.cpp @@ -444,284 +444,7 @@ computeBhpAtThpLimitProdWithAlq( DeferredLogger& deferred_logger, double alq_value) 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 = 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& 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& 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 range {controls.bhp_limit, *bhp_max}; - std::optional 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:: - 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 -std::optional -MultisegmentWellGeneric:: -bhpMax(const std::function& 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 -bool -MultisegmentWellGeneric:: -bisectBracket(const std::function& eq, - const std::array& range, - double& low, double& high, - std::optional& 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 -bool -MultisegmentWellGeneric:: -bruteForceBracket(const std::function& eq, - const std::array& 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; + return baseif_.computeBhpAtThpLimitProdCommon(frates, summary_state, maxPerfPress, rho, alq_value, deferred_logger); } template diff --git a/opm/simulators/wells/StandardWellGeneric.cpp b/opm/simulators/wells/StandardWellGeneric.cpp index c258c5814..53fa8d199 100644 --- a/opm/simulators/wells/StandardWellGeneric.cpp +++ b/opm/simulators/wells/StandardWellGeneric.cpp @@ -199,8 +199,13 @@ StandardWellGeneric:: computeBhpAtThpLimitProdWithAlq(const std::function(const double)>& frates, const SummaryState& summary_state, DeferredLogger& deferred_logger, + double maxPerfPress, 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 // rates and thp: // fbhp(rates, thp), diff --git a/opm/simulators/wells/StandardWellGeneric.hpp b/opm/simulators/wells/StandardWellGeneric.hpp index bdb18c297..1bdcb83bc 100644 --- a/opm/simulators/wells/StandardWellGeneric.hpp +++ b/opm/simulators/wells/StandardWellGeneric.hpp @@ -105,6 +105,7 @@ protected: std::optional computeBhpAtThpLimitProdWithAlq(const std::function(const double)>& frates, const SummaryState& summary_state, DeferredLogger& deferred_logger, + double maxPerfPress, double alq_value) const; // Base interface reference diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index c717c92ad..888257d71 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -2320,9 +2320,18 @@ namespace Opm 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::computeBhpAtThpLimitProdWithAlq(frates, summary_state, deferred_logger, + max_pressure, alq_value); } diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index aa1ffc772..dd40f439a 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -23,12 +23,14 @@ #include #include - +#include #include #include #include #include #include +#include +#include #include #include #include @@ -408,5 +410,288 @@ void WellInterfaceGeneric::reportWellSwitching(const SingleWellState& ws, Deferr deferred_logger.info(msg); } +std::optional +WellInterfaceGeneric:: +bhpMax(const std::function& 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& eq, + const std::array& range, + double& low, double& high, + std::optional& 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& eq, + const std::array& 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 +WellInterfaceGeneric:: +computeBhpAtThpLimitProdCommon(const std::function(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& 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& 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 range {controls.bhp_limit, *bhp_max}; + std::optional 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:: + 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 diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index b3666c5de..a775e8db0 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -177,6 +177,13 @@ public: bool changedToOpenThisStep() const { return this->changed_to_open_this_step_; } + std::optional computeBhpAtThpLimitProdCommon(const std::function(const double)>& frates, + const SummaryState& summary_state, + const double maxPerfPress, + const double rho, + const double alq_value, + DeferredLogger& deferred_logger + ) const; protected: bool getAllowCrossFlow() const; double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const; @@ -185,6 +192,23 @@ protected: WellTestState& well_test_state, DeferredLogger& deferred_logger) const; + std::optional bhpMax(const std::function& fflo, + const double bhp_limit, + const double maxPerfPress, + const double vfp_flo_front, + DeferredLogger& deferred_logger) const; + + bool bruteForceBracket(const std::function& eq, + const std::array& range, + double& low, double& high, + DeferredLogger& deferred_logger) const; + + bool bisectBracket(const std::function& eq, + const std::array& range, + double& low, double& high, + std::optional& approximate_solution, + DeferredLogger& deferred_logger) const; + // definition of the struct OperabilityStatus struct OperabilityStatus {