From 1001d35418e854ed14f778932c3865989945501c Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 9 Feb 2022 10:06:09 +0100 Subject: [PATCH 1/6] 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 { From ec08f80405d123f718dba07c4b92121bd18a45c8 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 9 Feb 2022 10:49:17 +0100 Subject: [PATCH 2/6] Refactor out the solving algorithm of bhp from thp via VFP to make it usable for injectors --- opm/simulators/wells/WellInterfaceGeneric.cpp | 28 ++++++++++++++++++- opm/simulators/wells/WellInterfaceGeneric.hpp | 10 +++++++ 2 files changed, 37 insertions(+), 1 deletion(-) diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index dd40f439a..a03b328d6 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -649,13 +649,39 @@ computeBhpAtThpLimitProdCommon(const std::function(const dou "find bhp-point where production becomes non-zero for well " + this->name()); return std::nullopt; } + const std::array range {controls.bhp_limit, *bhp_max}; + return computeBhpAtThpLimitCommon(frates, fbhp, range, deferred_logger); +} + +std::optional +WellInterfaceGeneric:: +computeBhpAtThpLimitCommon(const std::function(const double)>& frates, + const std::function)>& fbhp, + const std::array range, + DeferredLogger& deferred_logger) const +{ + // Given a VFP function returning bhp as a function of phase + // rates and thp: + // fbhp(rates, thp), + // a function extracting the particular flow rate used for VFP + // lookups: + // flo(rates) + // and the inflow function (assuming the reservoir is fixed): + // frates(bhp) + // we want to solve the equation: + // fbhp(frates(bhp, thplimit)) - bhp = 0 + // for bhp. + // + // This may result in 0, 1 or 2 solutions. If two solutions, + // the one corresponding to the lowest bhp (and therefore + // highest rate) should be returned. + // Define the equation we want to solve. auto eq = [&fbhp, &frates](double bhp) { return fbhp(frates(bhp)) - bhp; }; // Find appropriate brackets for the solution. - 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 diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index a775e8db0..53091a2a3 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -184,6 +184,9 @@ public: const double alq_value, DeferredLogger& deferred_logger ) const; + + + protected: bool getAllowCrossFlow() const; double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const; @@ -198,6 +201,13 @@ protected: const double vfp_flo_front, DeferredLogger& deferred_logger) const; + std::optional computeBhpAtThpLimitCommon( + const std::function(const double)>& frates, + const std::function)>& fbhp, + const std::array range, + DeferredLogger& deferred_logger) const; + + bool bruteForceBracket(const std::function& eq, const std::array& range, double& low, double& high, From e0573e99b11327febf87dcd62beb4b684d157e17 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 9 Feb 2022 10:59:24 +0100 Subject: [PATCH 3/6] iterate in computeBhpAtThpLimitProdWithAlq if no solution or potentials are negative --- opm/simulators/wells/StandardWell_impl.hpp | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 888257d71..f66e5b2a6 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -2328,11 +2328,30 @@ namespace Opm double pressure_cell = this->getPerfCellPressure(fs).value(); max_pressure = std::max(max_pressure, pressure_cell); } - return this->StandardWellGeneric::computeBhpAtThpLimitProdWithAlq(frates, + auto bhpAtLimit = this->StandardWellGeneric::computeBhpAtThpLimitProdWithAlq(frates, summary_state, deferred_logger, max_pressure, alq_value); + auto v = frates(*bhpAtLimit); + if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; })) + return bhpAtLimit; + + auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) { + // Solver the well iterations to see if we are + // able to get a solution with an update + // solution + std::vector rates(3); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger); + return rates; + }; + + return this->StandardWellGeneric::computeBhpAtThpLimitProdWithAlq(fratesIter, + summary_state, + deferred_logger, + max_pressure, + alq_value); + } From 650416c64759989d135de1e13c8639f05cb29947 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 10 Feb 2022 09:17:39 +0100 Subject: [PATCH 4/6] cleanup --- opm/simulators/wells/StandardWellGeneric.cpp | 189 +----------------- opm/simulators/wells/WellInterfaceGeneric.hpp | 4 +- 2 files changed, 3 insertions(+), 190 deletions(-) diff --git a/opm/simulators/wells/StandardWellGeneric.cpp b/opm/simulators/wells/StandardWellGeneric.cpp index 53fa8d199..e61618b4d 100644 --- a/opm/simulators/wells/StandardWellGeneric.cpp +++ b/opm/simulators/wells/StandardWellGeneric.cpp @@ -202,194 +202,7 @@ computeBhpAtThpLimitProdWithAlq(const std::function(const do 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), - // 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().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, getRho(), 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]); - }; - - // Get the flo samples, add extra samples at low rates and bhp - // limit point if necessary. Then the sign must be flipped - // since the VFP code expects that production flo values are - // negative. - std::vector 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); - } - for (double& x : flo_samples) { - x = -x; - } - - // Find bhp values for inflow relation corresponding to flo samples. - std::vector bhp_samples; - for (double flo_sample : flo_samples) { - if (flo_sample < -flo_bhp_limit) { - // We would have to go under 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 = 600.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() ? high : 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 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; - } + return baseif_.computeBhpAtThpLimitProdCommon(frates, summary_state, maxPerfPress, getRho(), alq_value, deferred_logger); } template diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 53091a2a3..e19811ae3 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -177,7 +177,7 @@ public: bool changedToOpenThisStep() const { return this->changed_to_open_this_step_; } - std::optional computeBhpAtThpLimitProdCommon(const std::function(const double)>& frates, + std::optional computeBhpAtThpLimitProdCommon(const std::function(const double)>& frates, const SummaryState& summary_state, const double maxPerfPress, const double rho, @@ -201,7 +201,7 @@ protected: const double vfp_flo_front, DeferredLogger& deferred_logger) const; - std::optional computeBhpAtThpLimitCommon( + std::optional computeBhpAtThpLimitCommon( const std::function(const double)>& frates, const std::function)>& fbhp, const std::array range, From cb99a2fc740d051d887a5a59e5a44c771aea9b29 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 11 Feb 2022 09:34:00 +0100 Subject: [PATCH 5/6] check also the validity of the iterated solution --- opm/simulators/wells/StandardWell_impl.hpp | 37 ++++++++++++---------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index f66e5b2a6..e0b1f8967 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -2333,25 +2333,30 @@ namespace Opm deferred_logger, max_pressure, alq_value); - auto v = frates(*bhpAtLimit); - if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; })) - return bhpAtLimit; + auto v = frates(*bhpAtLimit); + if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; })) + return bhpAtLimit; - auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) { - // Solver the well iterations to see if we are - // able to get a solution with an update - // solution - std::vector rates(3); - computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger); - return rates; - }; + auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) { + // Solver the well iterations to see if we are + // able to get a solution with an update + // solution + std::vector rates(3); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger); + return rates; + }; - return this->StandardWellGeneric::computeBhpAtThpLimitProdWithAlq(fratesIter, - summary_state, - deferred_logger, - max_pressure, - alq_value); + bhpAtLimit = this->StandardWellGeneric::computeBhpAtThpLimitProdWithAlq(fratesIter, + summary_state, + deferred_logger, + max_pressure, + alq_value); + v = frates(*bhpAtLimit); + if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; })) + return bhpAtLimit; + // we still don't get a valied solution. + return std::nullopt; } From 418880730ea6a301c8023060050c44ff73f240a9 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 15 Feb 2022 09:54:21 +0100 Subject: [PATCH 6/6] cleanup2 --- opm/simulators/wells/StandardWellGeneric.cpp | 2 +- opm/simulators/wells/WellInterfaceGeneric.cpp | 2 +- opm/simulators/wells/WellInterfaceGeneric.hpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/simulators/wells/StandardWellGeneric.cpp b/opm/simulators/wells/StandardWellGeneric.cpp index e61618b4d..a6de18d32 100644 --- a/opm/simulators/wells/StandardWellGeneric.cpp +++ b/opm/simulators/wells/StandardWellGeneric.cpp @@ -202,7 +202,7 @@ computeBhpAtThpLimitProdWithAlq(const std::function(const do double maxPerfPress, double alq_value) const { - return baseif_.computeBhpAtThpLimitProdCommon(frates, summary_state, maxPerfPress, getRho(), alq_value, deferred_logger); + return baseif_.computeBhpAtThpLimitProdCommon(frates, summary_state, maxPerfPress, getRho(), alq_value, deferred_logger); } template diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index a03b328d6..5e0d31506 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -657,7 +657,7 @@ std::optional WellInterfaceGeneric:: computeBhpAtThpLimitCommon(const std::function(const double)>& frates, const std::function)>& fbhp, - const std::array range, + const std::array& range, DeferredLogger& deferred_logger) const { // Given a VFP function returning bhp as a function of phase diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index e19811ae3..2459ff9ce 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -204,7 +204,7 @@ protected: std::optional computeBhpAtThpLimitCommon( const std::function(const double)>& frates, const std::function)>& fbhp, - const std::array range, + const std::array& range, DeferredLogger& deferred_logger) const;